FINDE: Neural Differential Equations for Finding and Preserving Invariant Quantities
Authors: Takashi Matsubara
Takaharu Yaguchi
Venue: ICLR 2023 poster (2023)
Keywords: neural ordinary differential equations
first integral
conservetaion law

Abstract: Many real-world dynamical systems are associated with first integrals (a.k.a. invariant quantities), which are quantities that remain unchanged over time. The discovery and understanding of first integrals are fundamental and important topics both in the natural sciences and in industrial applications. First integrals arise from the conservation laws of system energy, momentum, and mass, and from constraints on states; these are typically related to specific geometric structures of the governing equations. Existing neural networks designed to ensure such first integrals have shown excellent accuracy in modeling from data. However, these models incorporate the underlying structures, and in most situations where neural networks learn unknown systems, these structures are also unknown. This limitation needs to be overcome for scientific discovery and modeling of unknown systems. To this end, we propose first integral-preserving neural differential equation (FINDE). By leveraging the projection method and the discrete gradient method, FINDE finds and preserves first integrals from data, even in the absence of prior knowledge about underlying structures. Experimental results demonstrate that FINDE can predict future states of target systems much longer and find various quantities consistent with well-known first integrals in a unified manner.

Published as a conference paper at ICLR 2023
FINDE: NEURAL DIFFERENTIAL EQUATIONS FOR FINDING AND PRESERVING INVARIANT QUANTITIES
Takashi Matsubara
Takaharu Yaguchi
Osaka University
Kobe University
Toyonaka, Osaka, 560-8531 Japan
Kobe, Hyogo, 657-8501 Japan
mat matsubara@sss.eee es. www.ssu..kkaaaa.... acc. jp
yaguchi@pearl. kobe-u. ac. ip
ABSTRACT
Many real-world dynamical systems are associated with first integrals (a.k.a. in- variant quantities), which are quantities that remain unchanged over time. The discovery and understanding of first integrals are fundamental and important top- ics both in the natural seiences and in industrial applications. First integrals arise from the conservation laws of system energy, momentum, and mass, and from constraints on states; these are typically related to specific geometric structures the governing equations. Existing neural networks designed to ensure such first integrals have shown excellent accuracy in modeling from data. However, these models incorporate the underlying structures, and in most situations where neural networks learn unknown systems, these structures are also unknown. This limita- tion needs to be overcome for scientific discovery and modeling of unknown sys- tems. To this end, we propose firsi integral-preste preeervire neural differential equation (FINDE). By leveraging the projection method and the discrete gradient method, FINDE finds and preserves first integrals from data, even in the absence of prior knowledge about underlying structures. Experimental results demonstrate that FINDE can predict future states of target systems much longer and find various quantities consistent with well-known first integrals a unified manner.
INTRODUCTION
Modeling and predicting real-world systems are fundamental aspects of understanding the world in natural science and improving computer simulations in industry. Target systems include chemical dynamics for discovering new drugs (Raff et al.) (012)) climate dynamics for climate change predic-
tion and weather forecasting (Rasp et al. 2020; Trigo & Palutikof, 1999). and physical dynamics of
vehicles and robots for optimal control (Nelles 2001). addition to image processing and natural
language processing (Devlin et al. 2018 He et al. 2016). neural networks have been actively stud- ied for modeling dynamical systems (Nelles] 2001). Their history dates back to at least the 1990s
(see Chen al. (1990); Clouse et al. 1997); Levin & Narendral (1955); Narendra & Parthasarathy [9jb)e et al. (194)] Wang & Lin (998) for examples). Recently, two notable but dis- tinct families have been proposed. Physics-infone informeeddeee neural networks (PINNs) directly solve partial differential equations (PDEs) given as symbolic equations (Raissi et al. 2019). Neural ordinary dif- ferential equations (NODEs) learn ordinary differential equations (ODEs) from observed data and solve them using numerical integrators (Chen et al.) 2018). Our focus this time is on NODEs.
Most real-world systems are associated with first integrals (a.k.a. invariant quantities), which are quantities that remain unchanged over time (Hairer et al., 2006). First integrals arise from intrinsic geometric structures of systems and are sometimes more important than superficial dynamics in understanding systems (see Appendix[A for details). Many previous studies have extended NODEs by incorporating prior knowledge about first integrals and attempted to accurately learn a target system. Greydanus al. (2019) proposed the Hamiltonian neural network (HNN), which employs a neural network to approximate Hamilton's 's equation, thereby conserving the system energy called the Hamiltonian. Finzi et al. (2020a) proposed neural network architectures that conserve linear and angular momenta by utilizing the graph structure. Finzi et al. (2020b) also extended an HNN system with holonomic constraints, which led to first integrals such pendulum length.

---

Published as a conference paper at ICLR 2023
Table I: Comparison of Related Studies on Preservation of First Integrals.
Learning invariants Exact conservation
Energy
Monentum
Moss
Constraint
NODE (Chen al.
HNN Greydaauus
2019
LicConv Pinzi of ml 21p1m
DGNet Natubbra
CHNN
Tina et a A1400
NPM Yang et
Continuous HINDE (proposed)
Discrete FINDE (prppsed)
Matsubara et al. 2020) proposed a model that preserves the total mass of a discretized PDE. These studies have demonstrated that the more prior knowledge a neural network has about first integrals, the more accurate their dynamics prediction. See Table for comparisons.
Previous studies have mainly attempted to preserve known first integrals for better computer simula- tions. However, in situations where a neural network learns a target system, it is naturally expected that first integrals associated with the target system are unknown, and it is not clear which of the above methods are available. Therefore, this study proposes first integral-preservic preserving neural differ- ential equation (FINDE) to find and preserve unknown first integrals from data in unified manner. FINDE has two versions for continuous and discrete time; these have the following advantages.
Finding First Integrals Many studies have designed architectures or operations of neural net- works model continuous-time dynamics with known types of first integrals. However, the under- lying geometric structures of a target system are generally unknown in practice. In contrast, FINDE finds various types of first integrals from data in a unified manner and preserves them in predictions. For example, from an energy-dissipane FINDE can find not only known first integrals, but also unknown ones. Hence, FINDE can lead to scientific discoveries.
system, FINDE can find first integrals other than energy.
Combination with Known First Integrals FINDE can be combined with previously proposed neural networks designed to preserve known first integrals, such as HNNs. In addition, when some first integrals are known in advance, they can also be incorporated into FINDE to avoid rediscovery. Therefore, FINDE is available various situations.
Exact Preservation First Integrals The first integral associated with continuous-tin system is destroyed after the dynamics is temporally discretized for computer simulations. By leveraging the discrete gradient, the discrete-time version of FINDE preserves first integrals exactly (up to rounding errors) in discrete time and further improves the prediction performance.
2 BACKGROUND AND RELATED WORK
First Integrals Let us consider a time-invariant differential system u = f(u) on an N- dimensional manifold M, where denotes the system state and f M - T.M represents vector field on M. For simplicity, we suppose the manifold M to be a Fuclidean space RR. Definition 1 (first integral). A quantity V .M- + R is referred to d de u = f(u) if it remains constant along with any solution u(t). i.e. d
M
first integral of a system de V(u) = O.
4 a differential system de u = f(u) has K functionally independent first integrals Vi,..., Vk. the solution u(t) given initial value un stays at the (N K dimeeniinall submanifold
M' = fu M: Vi(u) = (uo)...... V*(u) = Vvx(ua))
(1)
The tangent space M' T of the submanifold M' M at point u is the orthogonal complement to the space spanned by the gradients VVa of the first integrals V₁ for k = 1,... K;
TuM- = {w E T.M: VVa(u) w 0 for k = 1,..., K}.
(2)
Conversely, the time-derivative o at point is on the tangent space T. for certain functions V 's the quantities VA 's are first integrals of the system u f(u); it holds that de Va (u)
=
(u) du VVe(u) f(u) 00
2

---

Published as a conference paper at ICLR 2023
One of the most well-known first integrals is the Hamiltonian H, which represents the system en- ergy of a Hamiltonian system. Noether's ' theorem states that a continuous symmetry of a system leads to a conservation law (and hence a first integral) (Hairer et al.] 2006). A Hamiltonian system is symmetric to translation in time, and the corresponding first integral is the Hamiltonian. Sym- metries to translation and rotation in space lead to the conservation of linear and angular momenta. However, not all first integrals are related to symmetries. A pendulum can be expressed in Cartesian coordinates, and then the rod length constrains the mass position. This type of constraint called a holonomic constraint and leads to first integrals. Models of disease spreads and chemical reactions have the total mass (population) as the first integral. Also for a system described by a PDE, the total mass is sometimes first integral (Furihata & Matsuo 2010). See Appendix A for the classes of dynamics, their geometric structures, and related studies to find or preserve first integrals.
First Integrals in Numerical Analysis For computer simulations, differential systems are dis- cretized in time and solved by numerical integration, causing numerical errors (which is composed of temporal discretization errors and rounding errors). Moreover, the geometric structures of the system are often destroyed, and the corresponding first integrals are no longer preserved. A com- mon remedy is symplectic integrator, which preserves the symplectic structure and accurately integrates Hamiltonian systems (Hairer et al., (006). However, the Marsden theorem states that symplectic integrator only approximately conserves the Hamiltonian (Zhong & Marsden 1988). Hence, many numerical schemes have also been investigated to preserve first integrals exactly, while these schemes cannot preserve the symplectic structure. Some examples are shown below.
Let the superscript E 1,..,..., f} denote the state u or time ts at the s-th time step, and AL" = = 1s denote time-step size. A projection method uses a numerical integrator predict the next state 7 from the current state u² and then projects the state ū onto the submanifold M' GGea 1986; Hairer et al., 2006 Section IV.4). The projected state *++ preserves the first integrals Vr. In particular, the projected state u*+1 is obtained by solving the optimization problem arg min u-+1
u-+1| - subject V& (u**))= Vu(") = 0 for k ,...,, K.
The local coordinate method defines a coordinate system on the neighborhood of the current state u² and integrates a differential equation on it (Potra & Yen| 1991 Hairer et al. 2006, Section IV:5). The discrete gradient method defines a discrete analogue to differential system and integrates it in discrete time, thereby preserving the Hamiltonian exactly (up to rounding errors) in discrete time (Furihata & Matsuo 2010: Gonzalez 1996 Hong et al. 2011).
Neural Networks to Preserve First Integrals NODE defines the right-hand side f of differen- f() using a neural network in the most general way with no associated first integrals (Chen et al., 2018.. NODE is universal approximator to ODEs and can approximate any ODE with arbitrary accuracy there an infinite amount of training data (Teshima et al., 2020). In practice, the amount of training data is limited, and prior knowledge about the target system is help- ful for learning (see [Sannai et al. (202T) for the case with convolutional neural networks (CNNs)). HNN (Greydanus et al.) 22019) assumes the target system to be a Hamiltonian system in the canon- ical form, thereby guarantecing various properties of Hamiltonian systems by definition, including the conservation of energy and preservation of the symplectic structure in continuous time (Hairer let al. 2006). Some studies have employed a symplectic integrator for HNN to preserve the energy and symplectic structure with smaller numerical errors (Chen al., 2020). LicConv and EMLP- HNN employ neural network architectures with translational and rotational symmetries to preserve momenta (Finzi et al. 2020a dynamics (Finzi et al., 2020b). Deep conservation extracts latent dynamics of a PDE system and preserves a quantity of interest by forcing its flux to be zero (Lee & Carlberg 202I). HNN++ also guarantees the conservation of mass in PDE systems by using coefficient matrix derived from dif- ferential operators (Matsubara et al. 2020). These methods preserve known types of first integrals and suffer from temporal discretization errors. In contrast, FINDE learns any types of first integrals from data and preserves them even after temporal discretization.
tial system 2 a
202T CHNN incorporates a known holonomic constraint in the
The neural projection method (NPM) learns fixed holonomic constraints using the projection (and inequality constraints) (Yang et al. 2020). DGNet employed discrete gradient methods to guar- antee the energy conservation in Hamiltonian systems (and the energy dissipation in friction sys- tems) (Matsubara et al.) 2020). While these methods preserve the aforementioned first integrals exactly in discrete time, their formulations are not available for other first integrals.
3

---

Published as a conference paper at ICLR 2023
Several studies have proposed neural networks to learn Lyapunov functions, which are expected to be -nncaance over time, contrast to first integrals (Manek & Kolter 2019: Takeishi & Kawahara 2020). If the state moves in the direction of increasing the function, it is projected onto or moved inside the contour line of the Lyapunov function. This concept is similar to that of the continuous-time version of FINDE but focuses
single non-increasing quantity in continuous time; FINDE preserves multiple quantities in both continuous and discrete time.
FIRST INTEGRAL-PRESEREE
NEURAL DIFFERENTIAL EQUATION
We suppose that a target system has least K unknown functionally independent first integrals. When a neural network learns the dynamics of the target system, it is not guaranteed to learn these first integrals. We suppose that certain neural network f for modeling the target dynamics is given, and in addition to this model f. we introduce a neural network that outputs K K-dimensional vector V(u) (V) V2(u)) VK ()) T. Each element is expected to learn one of the first integrals as V& R R for =1 ..... K. Then, the submanifold M' is defined in Eq. (1)
3.1 CONTINUOUS FINDE: TIME- DERIVATIVE PROJECTION METHOD
We propose a time-derivative projection method called continuous FINDE (cFINDE). The cFINDE projects the dimevverivtiive onto the tangent space M'.,' Roughly speaking, the cFINDE projects the dynamics on the space of the directions in which the first integrals do not change. In this way, the method can learn dynamics while preserving first integrals V, thereby finding unknown first integrals from data.
We refer to the neural network that defines the time-derivative f RRV - RR Applying the method of Lagrange multipliers to the projection method in Eq. (3), and taking the limit as the time-step size approaches zero, we have
RRN as the base model.
u f(u), f(u) f(u) M(u) A(u), at Vu))
where M = av and À E N is the Lagrange multiplier (see Appendix B.1 for detailed derivation). We transform the second equation obtain
By
d (u(())
av d du at u = ((w) M(w)(f() A(w)),
(5)
from which obtain the Lagrange multiplier A(u) (M(uM(w))))+M(()[ee( By elimi- nating A(u), we define the cFINDE as
du= f(u)
(())f((u) for (u) = M(M)MM((())MM (M) M(u).
(6)
Theorem (continuous-ting first integral preservation). The cFINDE 4 d
di u = f(u) preserves all dt Vv 0.
first integrals Vi k = 1,..,, K in continuous time, that is,
See Appendix B.1 for proof. The base model f can be NODE, HNN, any other model depending on the available prior knowledge. Additionally. if a first integral is already known, it can of the first integrals V₂ instead of being found by the neural network. Note
be directly used that even though the base model an HNN, due to the projection, the cFINDE f no longer Hamiltonian system in the strict sense.
Compared to the base model f. the cFINDE requires the additional computation of the neural net- work V. several matrix multiplications, and an inverse operation. The inverse operation has com- putational cost of O(K3), which is not costly if the number K of first integrals is small. Many previous models also need the inverse operation to satisfy the constraints and geometric structures, such Lagrangian neural network (LNN) (Cranmer et al. 2020). neural symplectic form (Chen ctt al. 2021). and CHNN (Finzi et al. 2020b).
3.2
DISCRETE FINDE: DISCRETE- TIME DERIVATIVE PROJECTION METHOD
The cFINDE is still ODE and hence needs to be solved using numerical integrator, which causes the temporal discretization errors in the first integrals. In order to eliminate these errors, is neces- sary constrain the destination (i.e., finite difference) rather than the direction (i.e., time-derivative).
4

---

Published as a conference paper at ICLR 2023
For this purpose, we propose discrete FINDE (dFINDE) by employing discrete gradients to define discrete tangent spaces, which are needed to constraint the state variables on the submanifold M
A discrete gradient VV is discrete analogue to gradient VV (Furihata & Matsuo 2010 Gonza- lez. 1996, Hong et al. 2011). Recall that a gradient VV of a function V RR + R can be regarded a function R
RV that satisfies the chain rule d V(u) VV(u) T d de
de u. Analogously
discrete gradient V is defined as follows:
Definition 2 (discrete gradient). discrete gradient VV of a function V RN - R is a function RN x - RN that satisfies V(v) V(u) VV (v, u) (u u) and V (u, VV(u).
The first condition is a discrete analogue to the chain rule when replacing the -ime-derivatives d V and with finite differences (V(v) - V(u)) and (v u). respectively, and the second condition ensures consistency with the ordinary gradient VV. A discrete gradient VV is not uniquely deter- mined and has been obtained manually. Recently. the automatic discrete differentiation algorithm (ADDA) has been proposed by Matsubara et al. (2020), which obtains a discrete gradient of a neural network in a manner similar to the automatic differentiation algorithm (Abadi et al. [0066] Paszke ct al. 2017). 2017). The discrete gradient is defined in discrete time; hence, the prediction using the discrete gradient is free from temporal discretization errors. See Appendix B.2 and the references Furihata & Matsuo 2010); Matsubara et al. (2020) for more details.
di
Following Christiansen et al. (2011); :): Dahlby et al. (2011), introduce discrete analogue to the tangent space Tur' called the discrete tangent space T(n.u)MM. In particular, for a pair of points (v, u) E MM,, the discrete tangent space is defined as T(w.u) {w E RR,, (V, u) T = for k I,...,K).
(7)
If the finite difference (u++1 u²) between the predicted and current states is on the discrete tangent space T( M',,,,, the first integrals V4 are preserved because Vx(u*+1)- Vu(u") =
(V+ =²" (u=*1 u²) 0. Note that similar concepts defined in different ways also
referred to as discrete tangent spaces (Cuell & Patrick 2009; Dehmamy et al. 2021).
We suppose that a neural network (c.g., NODE) f defines an ODE and numerical integrator pre- dicts the next state ã from a given state u². We call this process discrete-ting base model J. which satisfies --t = (ur,; At)). Subsequently, consider the model
AL"
xx ""11--" = (1u+11 ; A())
AP"
(ut uu; A1") ((uu; A1*) M(u*) uu"" A(wst) u²), V(u*-!) V(u") 0,
where M M(=+1 u²) (VVi(u*t) u²) VVk(u*t) w²)) As shown in Appendix B. this formulation is also derived from the projection method in Eq. Using the chain rule of the discrete gradient,
0 V(uett)- Mu u²) At
-'l-_ Muu+1 +1 u" (b u²; At"),
Ar"
Substituting this into Eq. (8) and eliminating the Lagrange multiplier X, define the dFINDE as -++ Ac" = ( uz; At*) = (I- (I-F(uu*h) ( u )))(") At") for Y Y = (MM M - M, (10)
where we have abbreviated M (u++1, u²)) and Y (wrt1, u²) to M and Y, respectively.
Theorem 2 (discrete-ting time first integral preservation). The dFINDE u"t1_u' = ( .....,K preserves all first integrals for k 1,..., K in discrete time, that is, Ve(u*+1)- Ve(ue) 0.
See Appendix B. for proof. Intuitively, dFINDE projects the finite difference (discrete-ting deriva- tive) et onto the discrete tangent space T after the numerical integration for each step, whereas cFINDE projects the time- -eme-derivative f onto the tangent space T,M' at every substep inside a numerical integrator. In the discrete-tima base model ve, the ODE f can be defined by any model, such as NODE or HNN, and the numerical integrator can be implemented by any method, such as the Runge-Kutte Kutta method or the leapfrog integrator. The projection method in Eq. (3). the method in Eq. (8). and the dFINDE in Eq (10 are implicit methods and hence relatively computationally expensive. However, only the dFINDE can be trained non- non-iteratively by standard backpropagation algorithms. As explained in Appendix B.3) this is because the next state u ++1 I is given during train- ing and the ADDA can explicitly obtain the discrete gradient and its computational graph.
5

---

Published as a conference paper at ICLR 2023
Table 2: Datasets, Dynamics, and First Integrals.
First Integrals
Dataset
Dynamics (Structure)
Energy Momentum
Mass
Constraint
Two-body problem
Canonical Hamiltonian
8
Discretized KdV oquation
Non-canonical Hamiltonian
50
Double pendulum
Poisson
8
FitzHugh-Nagumo model
Dirac
4 EXPERIMENTS
4.1 EXPERIMENTAL SETTINGS
Target Systems We evaluated FINDE and base models using datasets associated with first in- tegrals; these are summarized in Table 2 A gravitational two-body problem (2-body) on 2- dimensional configuration space is a typical Hamiltonian system in the canonical form. In addition to the total energy. the system has first integrals related to symmetries in space, namely, the linear and angular momenta. The Korteweg-De Vries (KdV) equation is a PDE model of shallow water waves. This equation is a Hamiltonian system in a non-canonical form and has the Hamiltonian, total mass, and many other quantities as first integrals. We discretized the KdV equation in space, obtaining ditteniiinii state u. A double pendulum (2-pend) is Hamiltonian system in po- lar coordinates. However, we transformed it to Cartesian coordinates; hence, it became a Poisson system. The lengths of the two rods work holonomic constraints and lead to four first integrals in addition to the Hamiltonian. The FitzHugh-Nagumo model is a biological neuron model as an electric eircuit, which exhibits a rapid and transient change of voltage called spike, As electric circuit, the currents through and voltages applied to the inductor and capacitor can be regarded as system states, which are constrained by the circuit topology and Kirchhoff's current and voltage laws. Then, this system has a state of four elements and two first integrals. Because the resistor dissipates the energy, the system not a Poisson system, but Dirac structure can be found (van) der Schaft & Jeltsemal 2014). We generated time-series set of each dataset with different initial conditions (hence, different values of first integrals). See Appendix foror more details.
Implementation We implemented the proposed FINDE and evaluated it under the following set- tings, We implemented all codes by modifying the officially released codes of HNN (Greydanus) al. 2019 and DGNet (Matsubara et al. 2020 We used Python 3.8. with packages scipy v. 1.7.3, pytorch 1.10.2, torchdiffeq v. 0. functorch v. 1.10 preview, and gplearn v. 0.4.2. We used the Dormand -ormand-Prince method (dopri5) (Dormand Prince 1986) as the numerical integrator, except in Section 4.2 All experiments were performed on single NVIDIA A100.
Following HNN Greydanus et al. 2019) and DGNet (Matsubara et al., 2020), we used fully- connected neural networks with two hidden layers. The input was the state u. and the output repre- sented the first integrals V for FINDE, time-derivative f for NODE, or the Hamiltonian H for HNN. Each hidden layer had 200 units and preceded a hyperbolic tangent activation function. Each weight matrix was initialized as an orthogonal matrix. For the KdV dataset, we used I-dimensional CNN, wherein the kernel size of each layer was 3. The double pendulum is a sccond-order system, imply- ing that the time-derivative az q of the position q is known to be the velocity Hence, we treated only the acceleration de v as the output to learn the 2-pend dataset. This assumption slightly improved the absolute performances but did not change the relative trends.
d
As the loss function for the cFINDE, we used the mean squared error (MSE) between the ground s+1 predicted from the current step uदr normalized truth future state uGT 1+1 and the future state u ppedd
by the time-step size At"; we named this the /-step error. For the dFINDE, we used the MSE between the left- and right-hand sides of Eq. 10p because the ground truth states u and *+1 are available during the training phase. The base model and FINDE were jointly trained using the Adam optimizer Kingma& & Ba 2015) with the parameters (31, 32) = (0.9, 0.999) and a batch size of 200.
htt httper///athhb.ooorrrrrre//////t

---

Published as a conference paper at ICLR 2023
The learning rate was initialized to 10 and decayed to zero with cosine annealing (Loshchilov & Hutter 2017) See Appendix B.3 and the enclosed source code for details about implementations.
Evaluation Metric We used the 1-step error as evaluation metric, which is identical the loss function for the cFINDE, and displayed it in the scale x x10-9. The lower this indicator, the better, as indicated by J. The MSES of the state system energy over long period are misleading indicators, as suggested in prior studies (Botev et al.] 2021; Jin et al. [020b; /laccaas Vlachas et al. (020). For example, a periodic orbit that is correctly learned except for a slight difference in angular velocity would have the same MSE an orbit that never moves from its initial position. Instead, we used the valid prediction time (VPT) (Botev et al. 2021 Jin al. 2020b) Vlachas al. (020). VPT denotes the time point 8 divided by the length 5 of the time-series at which the MSE of the predicted state u pred. first exceeds a given threshold 0 in an initial value problem, that is,
V PP (urped: MGT) = 3 max { [MSE (wered fred. "GT) < for all S5 05
(II)
The higher this indicator, the better, as indicated by 1. To obtain VPTs, we normalized each element of state to have the zero mean and unit variance in the training data and set 0 to 0.01. For systems with "spiking" behaviors, a small error in phase may be regarded as a significant error the state; for the Naguumo model, we obtained the VPTs by allowing for delay and advance of up
to steps.
analytical
leapfrog dFINDE
<FINDE
4.2 DEMONSTRATION OF FIRST INTEGRAL PRESERVATION
state
Before learning first integrals from data, we demonstrate that dFINDE can preserve first integrals without temporal discretiza- tion errors. We used a mass-spring system, which had the state u = (q v) dynamicss dr = v and at = -q, and system energy E(g,v) v) [0² + v²). Using an initial value of (1.0 0.0) and a time-step size of At = 0.2, we solved the initial value problem of the true ODE using the leapfrog integrator with or without FINDE, with the true system energy E as the first integral V. Notably, no neural networks nor training were involved.
0.0-
-1.0-
energy
0.500
d
0.495
step
Integration of a known mass-spring system by the leapfrog integrator. (top) States predicted by compari- son methods. (bottom) Energy calculated from the states pre- dicted.
Figure 1:
Figure shows the results, along with the analytical solution. The states predicted by comparison methods overlap and are apparently identical. However, the energy obtained by the leapfrog integrator fluctuates and the same is true for cFINDE. This is because the leapfrog integrator and cFINDE suffer from temporal discretization errors in first integrals. In con- trast, dFINDE preserves the energy accurately, the same the analytical solution. This is because dFINDE projects the state (q v) onto the discrete tangent space T(v.u).. M' at every step. Although smaller time-step size reduces temporal discretization errors, this result demonstrates the advantage of dFINDE. See Appendix D. for the case with the Dormand- Cormand-Prince integrator.
4.3 FINDING NON- NON-HAMILTONIAN FIRST INTEGRALS OF HAMILTONIAN SYSTEMS
We evaluated cFINDE and dFINDE on learning from the 2-body dataset. We used HNN as the base model f. We found that eFINDE and dFINDE obtained better performances if it did not treat the Hamiltonian H of the HNN one of the first integrals Vg. The medians and standard deviations of five trials are summarized in the leftmost column of Table 3 The cFINDE achieved better VPTs than the original HNN with K The dFINDE showed similar trend with slightly better performances; there is a trade-off between performance and computational cost. The HNN with either cFINDE or dFINDE found two first integrals in addition to the Hamiltonian H of the HNN. Even though a two-body problem is Hamiltonian system that an HNN can learn, the prior knowledge that there exist first integrals other than the Hamiltonian H can be clue that enables better learning. Despite their better long-term prediction performance, the HNN with either cFINDE or dFINDE yielded 1 1-step errors worse than the HNN, indicating that the 1-step error is misleading as evaluation criterion.
2, and its performance was suddenly degraded with K 3.
These example results are depicted in Fig. 2 (I), y1) and (x2, y2) became inaccurate in short time and the center-of-gravity position (Fe, yc)
In the absence of FINDE, the mass positions

---

Published as a conference paper at ICLR 2023
Table Results of cFINDE and dFINDE.
2-body HNN
Kdv
2-pend
FitzHugh-Nagumo
Model
K
1-step1
VPT1
I-step)
VPTT
1-stept
VPTI
1-stepl
VPTT
base model
5.17
0,362 20,026
5.59 20:30
0.339
=0038
0.82
2002
0. 10 =0035
73.66
0.236 20053
7.10
0.374
10.036
6.24:004:
0.371
-0088
0.75
:0:44
0.156
-0042
54.18
0.127
20
7.78
-1,39
0.450
=0.052
2.59
=0:11
0.608
-0,085
0.73
20.05
0.198
-0,008
37.03
0.437
00.044
+FFINDE
3
10
20.146"
3.19
=0:37
0.730
-0.991
0.69
=0,33
0.411
-0,03
>100
0.007 =0.007
4
10
0.101
20.005
3.65=0.30
0.641
=0071
0.77
=0.07
0.395
-0,83
5
10
0.080
20.014
4.68
20:43
0.601
-0.069
0.80
=0.07
0.585
-0.007
10
0.070
:0.019
7.79 =0:51
0.425
-0.067
12.53
20.00
0.005
-0.0000
1
7.01 1100
0.379
10,000
11.61
16.60
0,288
-0,083
0.75
10.10
0.152
00.017
47.07
18.03
0.117
10.122
2
7.03
0.475
=0,022
2.70
=0:26
0.598
=0,059
0.74
20.05
0.271
=0.111
33.24
23.40
0.455 00002
dFINDE
3
54.78 366.99
0.309
+0,024
3.78
+0.27
0.636
-0.024
0.69
#0.05
0.447
=0.081
319.70
491.11
0.049 +0.007
4
>> 10
0.102
t0.015
3.48
+0.32
0.780
-0,059
0.71
+0.03
0.454
-0.060
5
> 103
0.086
+0.011"
5.26
+0.15
0.718
-0.038
0.86
+0.09
0.591
-0.087
6
>> 103
0.059 20.017 9.60 #3.61 0.573 -0.121 58.88 222.98
0.037
=0,039
Standard deviation follows the + symbol; underlined results are better than those of the hase models; bold font indicates best results;" denotes trials that failed in training because of underflow of time-step size.
Notes:
true (stes Me))
[311: vn.)
22
(ste, Mo)
(122.22)
31.
Vec
0.5
2-
ground truth
HNN
+<FINDE
HNN
10a
+<FINDE 10
Figure/ 2: Example results of 2-body dataset. (left) Ground truth. (middle) HNN. (right) HNN with cFINDE.
Figure 3: Mean absolute errors for results of 2-body dataset. (left) HNN. (right) HNN with <FINDE.
21 Vi + ) deviated rapidly. The HNN with cFINDE accurately predicted the state for a longer period. Even after errors in the mass positions became non-negligible, errors in the center-of-gravity position were still small. Figure 3 shows the absolute errors averaged over all trials, which demon- strate how the trend changes with cFINDE. In both the x- and g-directions, the HNN without FINDE produced errors in the center-of oo-giiitty position ze (or yc), and those in the mass positions ø1. 22 (or y1, 322 at a similar level. In contrast, with the cFINDE, errors in the center-of- -ggavvittyy position were much smaller than those in the mass positions, implying that errors in one mass position can- celed out errors in the other. We performed a symbolic regression of first integrals V found by the neural network. For K = 2, the found first integrals V were identical to the linear momenta in the x- and g-directions up to affine transformation in most cases. See Appendix [.2 for detailed results. Therefore, we conclude that FINDE not only had better prediction accuracy but also found and preserved linear momenta (which are related to symmetries in space) more accurately despite not having prior knowledge about symmetries.
2
4.4 FINDING FIRST INTEGRALS OF UNKNOWN SYSTEMS
It is often unclear whether a target system is a Hamiltonian system or not, but one can expect that has several first integrals. We evaluated cFINDE and dFINDE using NODE as the base model and display the results in Table
For the KdV dataset, the NODE with either cFINDE or dFINDE obtained improved VPTs for a wide range of K. Agguwss an example result. The prediction states were apparently similar. In the absence of FINDE, the NODE increased all of its errors in proportion to time. With cFINDE, the error in total mass increased at the point where the two solitons collided, but then returned to the original level. Although the calculation is slightly inaccurate, the eFINDE learned to preserve the total mass. The error in energy continued to increase for K = 2, but remained within a small range for K 3. These results suggest that the first second quantity learned by the cFINDE was total
8

---

Published as a conference paper at ICLR 2023
50
50
50
Truth
NODE
+>FINDEE K=2
+<FINDE K=3
10
Truth
10-
0- NODE
0.5
K =
10
error in state
10 errorr in total mass 10*
error in energy 10
Figure 4: Example results of KdV dataset. (top) Predicted states. Red belts denote moving solitons.
(bottom) Mean absolute errors in states u, total mass k=1 Uk and energy, from left to right.
Truth
NODE
K=1
12
Ground Truth NODE +>FiNDÉ
Figure 5: Example results of 2-pend dataset for 2,000 steps. (left) Ground truth. (middle) NODE. (right) NODE
2000
2000
Figure 6: Example results of FitzHugh-Nagumo
with cFINDE for K 5.
dataset. Each panel shows one of four states.
mass, the third quantity was system energy, and the remaining quantity may correspond to one of the many first integrals of the KdV equation.
For the 2-pend dataset, the NODE with either cFINDE or dFINDE obtained improved VPTs with K 1 to 5. In addition to the system energy, the double pendulum has two holonomic constraints on the position, which lead to two additional constraints involving the velocity (see Appendix [ details). Thus, is reasonable that the NODE with either cFINDE or dFINDE obtained the best VPT for K = 5 first integrals and completely failed for K > first integrals. As exemplified in Fig- the NODE without FINDE did not preserve the lengths of rods, making the states deviate gradually. See Appendix D.3 for the case when actual constraints are known. For the FitzHugh-Naguma the NODE with either cFINDE or dFINDE obtained improved VPTs for K = 2. As exemplified in Fig. 6 the ground truth state converged periodic orbit, and only the NODE with cFINDE for K = reproduced similar dynamics. Without FINDE, the state did not remain in a limited region. wrong equilibrium; the sole quantity Vi may have attempted
dataset,
For K 1, the state converged and failed to learn both first integrals. We conclude that both cFINDE and dFINDE found all first integrals of the 2-pend and FitzHugh-Naguma datasets; K and K 2, respectively.
CONCLUSION
This study proposed first integral-preservnng neural differential equation (FINDE), which can find and preserve any type of first integrals from data in a unified manner. FINDE projects the time evolution onto the submanifold defined using the (discrete) gradients of first integrals represented by a neural network. We experimentally demonstrated that FINDE found and preserved first integrals that come from the energy and mass conservation laws, symmetries space, and constraints, thereby predicting the dynamics for far longer. FINDE is available even for an mnegy-ddssiipate system. When FINDE obtains the best prediction accuracy with K K², suggests that the target system has at least K' first integrals. Hence, FINDE has the potential to make scientific discoveries by revealing geometric structures of dynamical systems. See Appendix fo4 more discussions on K.
The numerical error tolerance 1 was negligible compared to the 1-step errors (which were 10-5 to 10 in absolute error). However, the dFINDE tended to obtain much better VPTs than the cFINDE. This result suggests that method leading smaller numerical errors produces model with smaller modeling errors, as observed in previous works (Chen et al., 2020 Matsubara et al, 2020). These results may form new frontier for integrating numerical and modeling emrors.

---

Published as a conference paper at ICLR 2023
REPRODUCIBILITY STATEMENT
See Section 4.1 for experimental settings. More detailed deseriptions can be found in Appendix B.3 for training procedure and Appendix [or datasets. The authors have enclosed the source code for generating the datasets and running the experiments as supplementary material.
ACKNOWLEDGEMENT
This study was partially supported by JST CREST (JPMICR1914), JST PRESTO (JPMIPR2IC7), and JSPS KAKENHI (19K20344, 20K11693).
REFERENCES
Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthicu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Rafal Jozefowic7, Yangqing Jia, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Mike Schuster, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large-scale machine learning on heterogeneous systems. USENIX Symposium Operating Systems Design and Implementation (OSDI), 2016.
Ferran Alet, Dylan Doblar, Allan Zhou, Joshua Tenenbaum, Kenji Kawaguchi, and Chelsea Finn. Noether Networks: Meta-Learning Useful Conserved Quantities, Advances Neural Information Processing Systems (NeurIPS), INeurIPS):1 1-00, 2021.
David G. T. Barrett and Benoit Dherin. Implicit Gradient Regularization. In International Confer- ence on Learning Representations (ICLR), 2021.
Aleksandar Botev, Andrew Jaegle, Peter Wirnsherger, Daniel Hennes, and Irina Higgins. Which priors materr? Benchmarking models for learning latent dynamics. In Advances in Neural Infor- mation Processing Systems (NeurIPS) Track on Datasets and Benchmarks, 2021.
Yuan Cao, Zhiying Fang. Yue Wu, Ding Xuan Zhou, and Quanquan Gu. Towards Understanding the Spectral Bias of Deep Learning. International Joint Conference on Artificial Intelligence (IJCAI), pp. 2205-2211, 2021.
Elena Celledoni, Andrea Leone. Davide Murari, and Brynjulf Owren. Learning Hamiltonians of Constrained Mechanical Systems. arXiv, pp. 1-18, 2022.
S. Chen, S. A. Billings, and P. M. Grant. Non-linear system identification using neural networks. International Journal of Control, 111611119111114, 1990.
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, David Duvenaud, Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems (NeurIPS), pp. 1-19, 2018.
Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi. Neural Symplectic Form Learning Hamil-
tonian Equations General Coordinate Systems. In Advances in Neural Information Processing Systems (NeurIPS), 2021.
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic Recurrent Neural Networks. International Conference on Learning Representations (ICLR), pp. 1-23, 2020.
Snorre Christiansen, Hans Munthe-Kass and Brynjulf Owren. Topics in strueture-preserving discretization. Acta Numerica, 20:1-119) 2011.
D S Clouse, C Giles, B G Home, and G W Cottrell. Time-delay neural networks: Representation and induction of finite-state machines. IEEE Transactions on Neural Networks, 8(5):1065-70, January 1997.
10

---

Published as a conference paper at ICLR 2023
Kevin L. Course, Trefor W. Evans, and Prasanth B. Nair. Weak form generalized Hamiltonian learning. In Advances in Neural Information Processing Systems (NeurIPS), number NeurIPS, 2020.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian Neural Networks. ICLR Deep Differential Equations Workshop, pp. 1-9, March 2020.
Charles Cuell and George W. Patrick. Geometric discrete analogues of tangent bundles and strained Lagrangian systems. Journal of Geometry and Physics, 59(7):976 997, 2009.
Morten Dahlby, Brynjulf Owren, and Takaharu Yaguchi. Preserving multiple first integrals by dis- crete gradients. Journal of Physics A: Mathematical and Theoretical, 44(30), 2011.
Nima Dehmamy, Robin Walters, Yanchen Liu, Dashun Wang, and Rose Yu. Automatic Symmetry Discovery with Lie Algebra Convolutional Network. In Advances in Neural Information Process- ing Systems (NeurIPS), number 2018, pp. 1-30, 2021.
Jacob Devlin, Ming- Wei Chang, Kenton Lee, and Kristina Toutanova. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. arXiv, pp. 1-15, October 2018.
R. Dormand and P. Prince. A reconsideration of some embedded Runge- Kutta formulae. Journal of Computational and Applied Mathematics, 15(2):203-211, 1986.
Mare Finzi, Samuel Stanton, Pavel Izmailov, and Andrew Gordon Wilson. Generalizing Convolu- tional Neural Networks for Equivariance to Lie Groups on Arbitrary Continuous Data. In Inter- national Conference Machine Learning (ICMI., pp. 3146-3157, 2020a.
Mare Finzi, Ke Alexander Wang, and Andrew Gordon Wilson. Simplifying Hamiltonian and La- grangian Neural Networks via Explicit Constraints. In Advances in Neural Information Process- ing Systems (NeurIPS), 2020b.
Mare Finzi, Max Welling, and Andrew Gordon Wilson. A Practical Method for Constructing Equiv- ariant Multilayer Perceptrons for Arbitrary Matrix Groups. In International Conference on Ma- chine Learning (ICML), 2021.
K. Fukunaga and D.R. Olsen. An Algorithm for Finding Intrinsic Dimensionality of Data. IEEE Transactions on Computers, 1-2012) February 1971.
Daisuke Furihata. A stable and conservative finite difference scheme for the Cahn- Hilliard equation. Numerische Mathematik, 87(4):675-699, February 2001.
Daisuke Furihata and Takayasu Matsuo. Discrete Variational Derivative Method: A Structure- Preserving Numerical Method for Partial Differential Equations. Chapman and Hall/CRC, De- cember 2010.
C. W. Gear. Maintaining Solution Invariants in the Numerical Solution of ODE SIAM Journal Scientific and Statistical Computing, 733):3447743 1986.
O. Gonzalez. Time integration and discrete Hamiltonian systems. Journal of Nonlinear Science, 46)4499466, September 1996.
Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian Neural Networks. In Advances in Neural Information Processing Systems (NeurIPS), pp. 1-16, 2019.
Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric Numerical Integration: Structure-
Preserving Algorithms for Ordinary Differential Equations, volume 31. Springer-Verre Berlin/Heidelberg, 2006.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition. In IEEE Conference Computer Vision and Pattern Recognition (CVPR), pp. 1-9, December 2016.
Philipp Holl, Nils Thuerey, and Vladlen Koltun. Learning to Control PDEs with Differentiable Physics. In International Conference Learning Representations (ICLR), 2020.

---

Published as a conference paper at ICLR 2023
Jialin Hong, Shuxing Zhai, and Jingjing Zhang. ferential Equations with a Conserved Quantity. 2017-2038, January 2011.
Discrete Gradient Approach to Stochastic Dif- SIAM Journal on Numerical Analysis, 49(5):
Eugene M. Izhikevich and Richard FitzHugh.
FitzHugh-Nagumo model.
http://scholarpedia.com htgg/aaticlo/iizHug---aaaee 2mogel, 2006.
Pengzhan Jin, Zhen Zhang, Ioannis G. Kevrekidis, and George Em Karniadakis. Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks. pp. 1-12, 2020a.
Pengzhan Jin, Aiqing Zhu, George Em Karniadakis, and Yifa Tang. Symplectic networks: Intrinsic prructure-rresevvee networks for identifying Hamiltonian systems. Neural Networks, 132:166- 179, 2020b.
Muhammad Firmansyah Kasim and Yi Heng Lim. Constants of motion network, August 2022.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In International Conference on Learning Representations (ICLR), pp. 1-15, December 2015.
Kookjin Lee and Kevin Carlberg. Deep Conservation: A latent-dynamics model for exact satisfac- tion of physical conservation laws. In AAAI Conference on Artificial Intelligence (AAAI), 2021.
Astiel U. Levin and Kumpati S. Narendra. Recursive identification using feedforward neural net- works. International Journal of Control, 6133:5333-547, 1995.
Ziming Liu and Max Tegmark. Machine Learning Conservation Laws from Trajectories. Physical Review Letters, 126(18):18066) May 2021.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-Net: Learning PDEs from Data. In International Conference on Machine Learning (ICML), pp. 3208-33033 3216. PMLR, July 2018.
Ilya Loshchilov and Frank Hutter. SGDR: Stochastic gradient descent with warm restarts. In Inter- national Conference on Learning Representations (ICLR), pp. 1-16, 2017.
Gaurav Manek and J. Zico Kolter. Learning Stable Deep Dynamics Models. In Advances in Neural Information Processing Systems (NeurlPS), pp. 1-9, 2019.
Takashi Matsubara, Ai Ishikawa, and Takaharu Yaguchi. Deep Energy-Based Modeling of Discrete- Time Physics. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
Robert M. Miura, Clifford S. Gardner, and Martin D. Kruskal. Korteweg-de Vries equation and gen- cralizations. II. Existence of conservation laws and constants of motion. Journal of Mathematical Physics, 18):12044 1968.
Kumpati S. Narendra and Kannan Parthasarathy. Identification and Control of Dynamical Systems Using Neural Networks. IEEE Transactions on Neural Networks, 1(1):4-27, 1990.
Oliver Nelles. Nonlinear System Identification. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001.
Adam Paszke, Gregory Chanan, Zeming Lin, Sam Gross, Edward Yang, Luca Antiga, and Zachary Devito. Automatic differentiation in PyTorch. Autodiff Workshop on Advances in Neural Information Processing Systems, pp. 1-4, 2017.
Florian A. Potra and Jeng Yen. Implicit numerical integration for euler-lagrange equations via tangent space parametrization. Mechanics of Structures and Machines, 19(1):77-98, 1991.
Lionel Raff, Ranga Komanduri, Martin Hagan, and Satish Bukkapatnam. Neural Networks in Chem- ical Reaction Dynamics. 2012.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informer neural networks: A deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
12

---

Published as a conference paper at ICLR 2023
Stephan Rasp. Peter D. Dueben, Sebastian Scher, Jonathan A. Weyn, Soukayna Mouatadid, and Nils Thuerey. WeatherBench: Benchmark Data Set for Data-Driven Weather Forecasting, Journal of Advances in Modeling Earth Systems, 12(11), 2020.
Akiyoshi Sannai, Masaaki Imaizumi, and Makoto Kawano. Improved Generalization Bounds of Group Invariant / Equivariant Deep Networks via Quotient Feature Spaces. In Conference on Uncertainty in Artificial Intelligence (UAI). October 2021.
Yujun Shen, Jinjin Gu, Xiaoou Tang, and Bolei Zhou. Interpreting the latent space of GANs for semantic face editing. Proceedings of the IEEE Computer Somputer Conference Computer Vision and Puttern Recognition, pp. 9240-9249, 2020.
Sjöberg, H. Hjalmarsson, and L. Ljung. Neural Networks in System Identification. IFAC Pro- ceedings Volumes, 27(8):359-33 382, 1994.
Yifan Sun, Linan Zhang, and Hayden Schaeffer. NeuPDE: Neural Network Based Ordinary and Partial Differential Equations for Modeling Time-Dependent Data. In Mathematical and Scientific Machine Learning Conference, pp. 352 372. PMLR, August 2020.
Naoya Takeishi and Yoshinobu Kawahara. Learning dynamics models with stable invariant sets. In
AAAI Conference Artificial Intelligence (AAAT), 2020.
Takeshi Teshima, Koichi Tojo, Masahiro Ikeda, Isao Ishikawa, and Kenta Oono. Universal Approxi- mation Property of Neural Ordinary Differential Equations. In NeurIPS Workshop on Differential Geometry Meets Deep Learning (DiffGeo4DI.)) 2020.
Ricardo M. Trigo and Jean P. Palutikof. Simulation of daily temperatures for climate change sce- narios over Portugal: A neural network model approach. Climate Research, 13(1):45-50 1999.
Arjan van der Schaft and Dimitri Jeltsema. Port-Hamiltonian Systems Theory: An Introductory Overview. Foundations and Trends@ in Systems and Control, 1(2):173-378, 2014.
R. Vlachas, J. Pathak, B. R. Hunt, T. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos. Backprop- agation algorithms and Reservoir Computing in Recurrent Neural Networks for the forecasting of complex spatiotemporal dynamics. Neural Networks, 126: 1611191217, 2020.
Yi Jen Wang and Chin Teng Lin. Runge-Kutte Kutta neural network for identification of dynamical sys- tems high accuracy. IEEE Transactions on Neural Networks, 9(2):294 307, 1998.
Shuqi Yang, Xingzhe He, and Bo Zhu. Learning Physical Constraints with Neural Projections. Advances in Neural Information Processing Systems, pp. 1-15, 2020.
Ge Zhong and Jerrold E Marsden. Lie-Poisson Hamilton-Jacobi Theory and Lie- iiePooisoon Integrators. Physics Letters A, 133(3):3-8, November 1988.
Yaofeng Desmond Zhong, Biswadip Dey. and Amit Chakraborty. Symplectic ODE-Net: Learning Hamiltonian Dynamics with Control. In International Conference on Learning Representations (ICLR), pp. 1 1-17, 2020a.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty, Dissipative SymODEN: Encoding Hamiltonian Dynamics with Dissipation and Control into Deep Learning. arXiv, pp. 1-6, 2020b.
13

---

Published as a conference paper at ICLR 2023
A HAMILTONIAN SYSTEM, ITS GENERALIZATION, AND FIRST INTEGRALS
Preliminary In this section, we briefly introduce potential target systems and related works. Meth- ods proposed by related works use specific prior knowledge about target systems, such constraints. In contrast, our proposed FINDE assumes a situation where neural networks learn systems with un- known properties. See, for example, Hairer et al. (2006); van der Schaft & Jeltsema (2014) for more details about geometric mechanics.
On an N N-dimensional manifold M, an ODE is defined using vector field f M - Tu which maps a point u on the manifold M to a tangent vector f(u) on the tangent space T..M. The NODE defines an ODE in this way (Chen et al. 2018). Given scalar-valued function H . - R on the manifold M, its differential d/f : M + TiM is a cotangent vector field (a.k.a. differential I-form), which maps point u on the manifold M cotangent vector d/ (u) on the cotangent space M.
Hamiltoanian System A Hamiltonian system is defined using a non-degenerate closed differential 2-form called symplectic form, which is a -ke-symmmetre bilinear map wa TuM TuM - at point u. A symplectic form assigned to manifold is called the symplectic structure. The coordinate-free form of Hamilton's 's equation is at for any w E Mu where XH is the Hamiltonian vector field. The symplectic form w gives rise to a bundle map wh T.MM + Mr,, with which Hamilton's is equation is rewritten as de u = (H())= = (dHH()). The right-hand side is locally equivalent to the product of coefficient matrix S and the gradient VH of the Hamiltonian H. Then, Hamilton's equation is obtained as de = SVH(u) Hamiltonian systems are often expressed in the canonical form, in other words, they are defined on Darboux coordinates, on which the state u the paired generalized position q and generalized momentum p. The corresponding coefficient matrix is S ( f for 2n N and the n- n-dimensional identity matrix L. The HNN was developed to model Hamiltonian systems in the canonical forms (Greydanus et al. 2019).
d
(u)),, wa wu(HH(u), w) (dll w)
d
An Euler-Lagrange equation with hyperregular Lagrangian and a Lotka-Voltera equation are also Hamiltonian systems; however, their coordinate systems are not Darboux coordinates. neural symplectic form (NSF) handles this class of equations (Chen et al. (02T). The KdV equation is also a Hamiltonian system not on Darboux coordinates. For Hamiltonian PDE systems, HNN++ was proposed (Matsubara et al., 2020). According to Darboux's 's theorem, any Hamiltonian system on an even-dim oven-dimensional manifold can be transformed into the canonical form.
Noether's theorem states that continuous symmetry of a system leads to conservation law. A Hamiltonian system symmetric (invariant) to translation in time and conserves the Hamiltonian H. A two-body problem is symmetric to translation and rotation in space and conserves linear and angular momenta. These quantities are first integrals. LieConv and EMLP-HNN had such symme- tries implemented in their architectures (Finzi al.) 2020a 202I). A pendulum is not symmetric to translation and rotation in space and does not conserve linear and angular momenta, but does exchange them with the base which it is fixed.
Poisson System A Poisson system is named after a Poisson bracket (. but it is convenient to refer it as a degenerate Hamiltonian system. A Poisson bracket is defined using a Poisson 2-vector B, which is a skmmesymeettic bilinear map Bu M x T - R at point u. The Poisson 2-vector gives rise to bundle map B T M + TrM and defines Hamilton's equation as de u = (u)H. The Darboux-Lie theorem states that any Poisson system can be transformed
d
into the canonical form fu = SVH(w) by using a matrix 5 = (-r 0 for 2h < N. The last
dt 0
0 0
N 2k elements remain unchanged and correspond to the first integrals. In this sense, a Poisson system is degenerate Hamiltonian system. A Poisson 2-vector assigned to a manifold is called Poisson structure. Several models of the dynamics of disease spreading and chemical reactions are Poisson systems, and total population and molecular mass are typical first integrals.
A Poisson neural network (PNN) learns to transform a given Poisson system into a canonical form (Jin al. 2020a).
14

---

Published as a conference paper at ICLR 2023
Constrained Hamiltonian System A constraint C(q) = 0 on the position q is called holonomic constraint. Holonomic constraint appear, for example, when the arm's 's length restricts the position of robot's 's hand. Differentiating holonomic constraint C(q) 0 yields constraint involving the velocity G(g, v) 8g constraint leads to two first integrals C' and G. A Hamiltonian system with holonomic constraints is also a Poisson system; in particular, it is a constrained Hamiltonian system.
ac O, which is simply called velocity constraint. Hence, each holonomic
A CHNN incorporates the known holonomic constraints C(q) and corresponding velocity con- straints Gig, v) of a Hamiltonian system in the canonical form (Finzi et al. 2020b). The origi- nal study suggested that CHNN may learn holonomic constraints from data, but this has not been tested. For modeling a constrained Hamiltonian system, it is sufficient to incorporate only velocity constraints G(q. a) because a holonomic constraint C(q) is implicitly satisfied if the corresponding velocity constraint G(q. v) is satisfied. Celledoni et al.) (0022) used such formulation, and extended HNN and CHNN to systems on non-Euclidean spaces. A neural projection method learns fixed holo- nomic constraints, as well as inequality constraints, which are outside the scope of this study (Yang) al., 2020). This method updates the state by solving an optimization problem similar to Eq. (3) iteratively using the gradient descent method at every training step. Subsequently, it applies the backpropagation algorithm to all the optimization iterations. Thus, it has high computational and memory costs.
These studies mainly focused on sisiicaall-inducce holonomic constraints and may not work for other first integrals, shown in Appendices[D.] D.5] However, the purpose of FINDE is to find and preserve general first integrals, including energy and mass not limited by constraints.
Dirac Structure A Dirac structure is named after Dirac bracket, generalization of the Poisson bracket (van der Schaft & Jeltsema 2014), and can be found in various systems. For a rolling disk, the direction in which the disk can move forward without slipping is limited by the disk's orien- tation. This constraint is called a non-holonomic constraint. In an electric circuit, when elements are connected in series, the current flow through each element always the same. This constraint is called Kirchhoff's s current law. One can find Dirac structures in these systems. The dissipative SymODEN was proposed to model a port- Horr-tmmltonn system the canonical form (Zhong al.) 2020b), which is a special case of the Dirac structure. To the best of our knowledge, a neural network model for a general Dirac structure has not yet been proposed. FINDE is the first neural network method to learn Dirac structures better than NODE can, even though it is not specialized for Dirac structures.
PDE with Mass Conservation The total mass of a PDE system is sometimes preserved (Furihata & Matsuo 2010). The KdV equation is a Hamiltonian system that describes shallow water waves, in which the energy and total mass preserved. The Cahn-Hilliard equation is a model of phase separation of copolymer melts, in which the total mass is preserved, but the energy is dissipated. In general, quantity in an area is preserved if its flux entering minus its flux leaving is zero. Deep conservation extracts latent dynamics of a PDE system and preserves a quantity of interest by forcing its flux be zero (Lee & Carlberg 2021). ). HNN++ also ensures mass conservation by designing coefficient matrix that determines local interaction (Matsubara et al. 2020).
General First Intergals A concurrent study, "Constants-of-moote penalty loss function so that NODES learn to preserve first integrals (Kasim & Lim 2022); how- ever, unlike other related methods, this method does not guarantee preservation. A Noether network was proposed to model videos that do not always capture physical phenomena (Alet et al.) 2021). A subset of the latent variable is assumed to represent image features that do not change during a video, such as the appearance of objects. For prediction, these features are forced not to change. The Noether network potentially useful for learning physical phenomena from videos, but is more similar to semantic manipulation of latent variables (Shen et al. 2020).
network," introduced the
Some studies have investigated methods that do not predict dynamics but specialize in finding first integrals (Fukunaga & Olsen, 1971 Liu & Tegmark 2021). These methods can be used to help FINDE determine the hyperparameter K. They commonly estimate the number (N- - K) of dimen- sions of the tangent space M'' of the submanifold M' at point u using its neighbors. For example, AI Poincaré proposed by [Lu & Tegmark (202T) assumes that all data points share the submanifold M' and uses autoencoder to reconstruct the tangent space T..M.. Hence, it can only process
15

---

Published as a conference paper at ICLR 2023
single long time series with fixed first integrals. In contrast, our proposed FINDE can leverage dataset multiple time series with different values the first integrals.
B DETAILS OF METHODS
B.1 DERIVATION OF FINDE
Continuous FINDE (cFINDE) Let u denote current state and f denote vector field. After a time interval At the state transitions to u++1 typical projection method projects the state a+ onto a submanifold M' and obtains state u=-1, which preserves the first integrals V
(V1 Vx) This procedure is defined as optimization problem in Eq. (
arg min 2 u subject to Vi(u*t1) Vi(u") = 0 for k = L....,K K.
(AI)
One can solve the problem using the method of Lagrange multipliers. A Lagrangian function is F(u- X) || 1=11 + (V(u*+1) V(u))
(A2)
where X is the Lagrange multiplier. The stationary point satisfies
OF
T u ut + Tu
X = 0,
Bu-
(A3)
OF V(u++)) Vuu)) =.
8x
Subsequently, projection method can be redefined as
u s+1 ū -++1
av y X, (u
(A4)
V(u*+1) V(u") 0.
We transform Eq. A4)
into
11
11 n²"
av
X,
At
Ae
(A5)
V(u +1 V(u")
0,
A1
where À X/At. Taking the limit At > +0, we obtain Eq.
f(u") f(u²) (owr A at (u")
(A6)
The second equation ensures that a state transition following the new vector field preserves the first integrals V. By climinating the Lagrange multiplier A(u), we define the cFINDE as in Eq. (6). that is,
de u f(u) = (u))( for Y(u) = M(u) ((M)T((M())) (M)
(A7)
where M av version of projection method. The preservation of first integrals can be proved follows.
du. Because of the above derivation, the cFINDE can be considered a continuous-time
Proof of Theorem
7
d (u)
V V(u)
aV u du di
dt
(())((u)
MTu)(( MM)/ MMM (M(u) - (())M (M(u) - M(w))f(u)
0.
Hence, holds that 4 V%V(w) (w) = for = 1.... K, indicating that the cFINDE du u = f(u) preserves all first integrals Ve in continuous time.
16

---

Published as a conference paper at ICLR 2023
Discrete FINDE (dFINDE) For dFINDE, we take the discrete gradient of the Lagrangian equa- tion in Eq. A2 and obtain the discrete version of the necessary conditions for the stationary point;
(u Luy = u ū+ + Mu +1 u ) = O,
(A8)
OF ox
= v(uut) V(u") (u")== -.
M(u+1 u²) corresponds to the Jacobian <uu; At") and the dFINDE u*+1
#V Zu
By substituting the base model 7
At"
(u u² A(*) into the above equation and dividing
Are
the first equation by At, we obtain Eq. 8
n+
A0 v uu; At"),
(u(ut1 u² A1)) = A(") - M(u ([u+tt,"))) uu) A(u ,u),,
(A9)
V(u*t)- V(u")= =0
where À= 1// t. By eliminating the Lagrange multiplier X, we define the dFINDE as in Eq.
10
that is,
1 -² (u u²; At") (1 Yu u )),(" At") for Y M (M M ) )
Ar
(A10)
The preservation of first integrals can be proved as follows.
Proof of Theorem 2
V(u V(u') = M Mu++1 u² uu) u²)
M(u+1 ut" (bb( u²; AF)
M( - MI' (M M ) M ( u² At
(M (MM )( M Ma Att At
(M M M) +1 u"; At
0.
Hence, it holds that V/ V2(=*t) Vi(u") for k: = 1,..., ., K. indicating that the dFINDE (b(++1 u²; At") preserves all first integrals V& in discrete time.
u
B.2 DISCRETE GRADIENT
discrete gradient is discrete analogue to gradient (Furihata Matsuo 2010 Gonzalez 1996
Hong et al. 2011). Discrete gradients that satisfy Deñnition 2 are not unique, and many varia- tions have been proposed. For neural network, Matsubara et al. (2020) proposed the automatic discrete differentiation algorithm (ADDA), We briefly introduce the algorithm in the case of finite- dimensional Euclidean spaces. The differential dg of function RA R&M a linear operator dgu: RN RM at point and satisfies
llg(u + h)) g(u)= + dgu (h))ll 0. |
lim
(A11)
hh |Rx
The differential dg acting on a vector w is equivalent to the product of a vector w with the Jacobian Jy(u) of the function g at point u: dgu (w) Jatn) w. Similarly, according to the chain rule, the differential d(h o g) of composition h g of functions g, h is equivalent the multiplication with a series Ja(g(u)) ggw) of Jacobians. Therefore, the automatic differentiation algorithm obtains the differential of a neural network. The differential dg of a function g N + IR is horizontal vector, and the gradient Vg of the function is a vertical vector dual to the differential. Therefore, the gradient Vg is obtained by transposing the differential dg. The ADDA replaces each Jacobian with its discrete analogue. For linear layers, such as fully-connected and convolution layers, the discrete Jacobian is identical to the ordinary Jacobian. For element-win nonlinear layers, such as activation functions, a diagonal matrix composed the slopes between two inputs can act as the discrete Jacobian. A discrete gradient obtained by the above steps satisfies Definition] 2
17

---

Published as a conference paper at ICLR 2023
B.3 PREDICTION AND TRAINING PROCEDURES
For ODEs modeled by neural networks. various training and prediction strategies have been pro- posed to date (Chen et al.] 2018 2020: Course et al. 2020 Matsubara et al. 2020) Zhong et al.) 2020a); FINDE can adopt any of these. In our experiments, used the following simple strategies.
In the case of the cFINDE and base models, taking a state &ar from the dataset, a numerical inte- grator solves the ODE
s+1 This process can be informally d u = f(v) and predicts the next state pred.
expressed as at
-1- /
u
s+1 pred. = u +
[(u((T)dd for u((").
(A 12)
We solved this integration using prchhiifeq.odeimg The prediction accuracy can be evaluated using the difference between the predicted state upred and ground truth ust taken from the dataset. We normalized the difference by the time-step size At and defined the I-step error Lisstop as
++1 sit Ci-swp( ( uped UGT
uGT. At")
s+1 UGT
u
u 6+1
pred.
u
(A13)
At"
At"
The cFINDE and base models were trained to minimize the 1-step error Li-step-
In the case of the dFINDE, the next state "pedd is predicted by solving Eq. (10 an implicit scheme; in particular,
s+1 pred.
u
arg min
s+1 (uupti u (Z Atr)
(A14)
At
Therefore, prediction by the dFINDE is implicit. For evaluation, we solved this scheme using scipy. pptime and obtained the 1-step error in Eq. (A13). However, during the training phase, the ground truth uar of the next state known. Hence, we substituted this into Eq. (10), and then used the difference between the left- and right-hand sides of the dFINDE the loss function:
Cu #+1 UGT, Atr)
21
s+1 GT
u
(
s+1 uGT)) A
(A15)
A
The discrete Jacobian M (and hence Y) can be obtained explicitly, and an explicit numerical inte- grator can be used for the base model vb. Hence, the process to obtain the value of the loss function explicit, and the dFINDE can be trained an explicit way, whereas the prediction is still implicit.
Some previous studies have proposed alternative strategies. For example, a loss function can be defined as the sum of the errors at multiple time points during long-term prediction. The cFINDE naturally adopt such a training strategy, and the dFINDE can adopt after minor modification. While it is helpful to pursue absolute performance, requires additional hyperparameters, such as the length of prediction time, and additional effort adjust them. We used the 1-step error in the present study for simplicity and fair comparisons.
The function V(u) learning first integral may become constant function during training; subse- quently, its Jacobian matrix vanishes (v = 0). this case, our algorithm returns a division-by-
3u
zero error because requires the inverse of the matrix (V( av(u) T for the projection. We have not taken any special measures to prevent such errors, but no errors occurred in any experiments with proper settings. The division-by-zen
Du du
eπrоrs have occurred only when FINDE assumes unreason- able number of first integrals (c.g., K for the double pendulum, which has five first integrals). FINDE works correctly even when the functions f(u) and V(u) learn the same first integrals; verified such a case Section 4.2 where both functions are known.
FINDE learns first integrals point-by-point and the found first integral is not always consistent over the domain. The same can said about the energy function of HNN, and this type of problem an open problem for neural network models of dynamical systems.
18

---

Published as a conference paper at ICLR 2023
C DETAILS OF DATASETS
To generate each dataset, we used scipy package and the Normand-Prinet method (dopri5) with the default relative tolerance of 10-9, unless otherwise stated. Experiments on the KdV dataset were performed with double precision, and all other experiments were performed with single precision.
Hamiltonian System in Canonical Form: Two-Body Problem A gravitational two-body prob- lem on a 2-dimensional configuration space has state u composed of the 4 4-dimensional position q = (Ø1 V1 22 3/2 and 4-dimensional velocity v = (va1 Vy1 U22 Uy2)". T. This is a second- di 4 = v. The momentum Px, I1 equals m1 The time- dt of the velocity v is called the acceleration. The acceleration of I is given by de Uz, -mmm
order ODE, indicating that d
derivative d d
(( ((=1-22)² + -2 32) )ay3/2, where G, m,, and m2 denote the constant of gravity
-1-22
and masses of two bodies, respectively. The same process applies for the remaining positions.
The total energy of the two-body problem is given by
Gm1m²
H = - [mi (c + v + +2 + x2)) V(n (ic1-33) x2) (y) y2)²
(A16)
The first and second terms denote the kinetic and potential energies, respectively. The two-body problem is Hamiltonian system, and the dynamics mentioned above can be rewritten as Hamilton's equation. The Hamiltonian H is one of the first integrals; the two-body problem has other first integrals, such as the linear momenta in the x- and y-directions
M1 + M2022
M1 + 122
Ps
Py
(A17)
my + m2
mi m2
and angular momentum (Hairer et 2006).
We set G, mi, and m2 to 1.0. The initial distance Tı V of + of of mass my from the origin was set to r1 (0.5,, 1.0), and the initial angle 01 tan -()) was set to 01 4(0, 2T). The
initial speed |/ = / 02 + u was set to 2 Cu, where G, ~ (1,, 0.05). The initial angle of the velocity was set to 0 + 0.5x + EOT, where Eo ~ N' 0.05). The initial condition of the other mass my was set to the opposite of the mass m1. Subsequently. the two masses trace elliptical orbits, and when E = Eg = 0. they trace exactly circular orbits. In addition, we added a perturbation following (O, 0.01) to the velocities of both masses, which corresponds to the center-of -gavitty velocity.
We set the time-step size At to 0.01 and generated 1,000 time-series of 5 = 500 steps for training and 10 time-series of 5 10, 000 steps for evaluation. We trained each model for 100,000 iterations.
Hamiltonian System in Non-Canonical Form: KdV equation The KdV equation is model of shallow water waves and is known have soliton solutions (Furihata) 2001). The dynamics given by
Mt -uul + Burrz,
(A18)
where x denotes the spatial position and the subscripts denote partial derivatives; for example, ut 8u The Hamiltonian is given by
1 au H(u)= -
1 2 33² dr.
(A19)
As Hamilton's equation d a SVH, the partial differential operator dz acts as the coefficient matrix S. This system is [iouville integrable and has infinitely many first integrals, including the Hamiltonian H, total mass 4₁ f udæ, and T2 f u²dx (Miura et al. 1968 Other first integrals are defined using higher-order partial derivatives.
For PDEs, PINNs are known to provide solutions when symbolic equations and boundary conditions are given (Raissi et al.) 2019). We, in contrast, consider learning spatially discretized PDEs as ODEs from observed data and solving them using numerical integrators, in the same context as NODEs and HNNs; this topic has also been studied extensively (Long et al, 2018) Matsubara et al., 2020 Holl et al. 2020). Following the experiments in a previous study (Matsubara al.
Sun al. 2020
19

---

Published as a conference paper at ICLR 2023
2020), we discretized the KdV equation in space; it no longer has infinitely many first integrals. We set a = -6, 3 = 1, spatial size to 10 space units, and space mesh size to 0.2: the system state had 50 elements. We generated two solitons as the initial condition; each was expressed as 12 x (x(x d)), where the size followed (0,55, 2) and the initial position d of one soliton was set to be at least 2.0 from that the other.
We set the time-step size At to 0.001 and generated 1,000 time-series of S 500 steps for training and 10 time-series of = 10, 000 steps for evaluation, using the discrete gradient method to ensure energy conservation (Furihata 2001). We trained each model for 30,000 iterations.
Due to the spatial discretization, the KdV dataset contains spatial truncation errors. When the neural network learns this dataset, spatial truncation errors are additionally introduced. An evaluation using the analytical solution as a dataset or datasets created with different spatial resolutions is included in future work.
Poisson System: Double Pendulum A double pendulum (2-pend) is depicted in Fig. Al In polar coordinates, this is a Hamiltonian system. The state is composed of the angles (01. U2) of the two rods and their angular velocities (w1, w2). This is also second-order ODE, indicating that de
d 01 and d di 82 = w2.
Let 4, la denote the lengths of the two rods, m1. m2 denote the masses of the two weights, and g denote the gravitational acceleration. The acceleration is given by
d W1 m2g sin O2 cos A (11wf cos A + (2w2 mm sin A (m² m2)g sin 0g
di
(1mm + m2 sin² A)
(A20)
dt W2 = (m + (huw sin A gsin 02 + y sin 02 cos A) + m2l2w3 sin cos
12(m1 m2 sin? A)
where = 02. In 2-dimensional Cartesian coordinates, the state is composed of the positions (I1, y1, 22 22, V2) of the two masses and the corresponding velocities (Uz1, Uyl, U22, Uy2). The position is transformed by #1 = 4 sin 01. f/1 = l cos 01. 22 = Ii + l2 sin 02. and y/2 = 8/1 - l2 cos 02. and the velocity is transformed accordingly. The total energy H is given by
11 - ( - " m2(e2 )) )) + g(miv) m232).
(A21)
The first and second terms denote the kinetic and potential energies, respec- tively. The double pendulum is no longer Hamiltonian system in Carte- sian coordinates. Because the lengths of the two rods are constant, the double pendulum has two constraints on the position: If = of + vी and i2 = (I2 - 2 + (V2 ) These constraints are holonomic constraints, and they lead to constraints involving the velocity, namely 0 = I1 U2, + 1/ and 0 = (x22 (tx, Uz, (3/2 y1 )1) v. ). When the constraints involving the velocity are satisfied, the holonomic constraints are implicitly satisfied. Therefore, the number of first integrals is five; however, three first integrals are sufficient to determine the dynamics. The dynamics degener- and classified as a constrained Hamiltonian system, or a Poisson system in more general case.
y
x
A
Ug
l2
m
0z
m2
Figure A1:
Dia-
gram of the double pendulum.
We the masses of the weights mi m2 1.0 and the gravitational acceleration g to 9.8. We set the lengths 4, l2 of the two rods to follow (0.9, 1.1), the initial angles 0,, O2 to follow u(-0.5, 0.5), and the initial angular velocities 61, 8g to follow ((.10.1,0.))
We the time-step size At to 0.1 and generated 1,000 time-series of 500 steps for training and 10 time-series of S 5, 000 steps for evaluation. We trained each model for 100,000 iterations.
Dirac Structure: -Naguuo Model FitzHugh proposed a model of the electrical dynamics of a biological neuron, and J. Nagumo created equivalent electric circuit. This model called the the van der Pol oscillator, the state oscillates when the magnitude of the external current source within an appropriate range. The circuit comprises a resistor R, inductor L, capacitor C, tunnel diode D, and voltage source E connected shown in Fig. A2 The whole circuit is connected to
model (Izhikevich & FitzHugh) 2006) and is a modified version of
20

---

Published as a conference paper at ICLR 2023
an external current source I. Let In denote the current through the resistor R. and VR denote the applied voltage. Ohm's 's law and other properties of the elements give VR = IRR, C de Vc = Ics L J = Vz. and TD = D(VD), where we treat D as a nonlinear function. Kirchhoff's current law (KCL) obtains Ic ID + In I and In I, and Kirchhoff's voltage law (KVL) obtains I and V = Ve. and set L 1/0.08, R = 0.8, C. 1.0, VE -0.7, and D(V) = V3/3 V. Subsequently, we obtain the FiizHugh-Nagumo
Vc - VD VR Vt, + E. We denote W
model of the original parameters
d
V V 13/3 W +1,
dt
(A22)
d
W 0.08(V 0.7 0.8W).
dt
Due to the resistor R. the BertHugg-Nagunne model is not energy-conserving system.
Consider a situation where the current through and the voltage ap- plied to stateful elements (capacitors and inductors) are measurable, but the connections between the elements are unknown. We treated Ic, Ir. Vc, VT. as the system state u. Because the state is in 4- dimensional space and the dynamics is intrinsically 2-dimensional, there exist two first integrals; for example, but not limited to. I To + D(Vc) + I and = Ve J Vz. This type of elec- tric circuit is example of Dirac structure because the state vari- ables are constrained by the circuit topology and Kirchhoff's s cur- rent and voltage laws (van der Schaft & Jeltsema, 2014). From the viewpoint of generalized Hamiltonian systems, (Z,, Vc) corre- sponds the position, and (VL, Ic) corresponds the momentum. The electric circuit can be described as a nort-Hamiltonian system in a non-canonical form. Because of the non-canonical form, the FizzHug model is outside the scope of CHNN and dissi- pative SymODEN (Finzi et al. 2020b Zhong et al. 2020b).
R
L
E
Figure A2: Circuit dia- gram of FitzHugh-Nagumo model Izhikevich & FitzHugh 2006).
We set the external current source / to follow (0.77, 1.1), set the initial values of V and W to follow (1-1.5, 1.5) and (0.0, 2.0), and transformed them to the state.
We set the time-step size A to 0.1 and generated 1,000 time-series of S = 500 steps for training and time-series of 2, 000 steps for evaluation. We trained each model for 30,000 iterations.
D ADDITIONAL RESULTS AND DISCUSSION
analytical cFINDE
dopri5
D.1 DEMONSTRATION OF FIRST INTEGRAL PRESERVATION
dFINDE
energy state MSE
In Fig. we examined a mass-spring system and FINDE using the leapfrog integrator. We also examined the case with the Dormand- Prince integrator (dopri5), as shown in Fig. A3 We increased the
10
0-
number of steps to 10a, and displayed the MSEs of the state in- stead of the state itself. First, we focus on the energy. Even using the Dormand-Prince integrator, fourth-order method, the energy is slightly decreased. The cFINDE with the Hormand-Prinee inte- grator shows the same tendency. This phenomenon is due tem- poral discretization errors and is called energy drift. The dFINDE with the Dormand-Prien Prince integrator significantly suppresses the error in energy. The remaining error is caused by rounding errors.
0.500
0.499
step
10a
Figure A3: Integration of known mass-spring system by Normand-Prine integrator. (top) Mean squared errors in states predicted by compari- son methods. (bottom) Energy calculated from the states pre- dicted.
When the focus is on the MSEs of the state, the trend is different: the dFINDE with the Dormand-Prince integrator suffers from the most significant errors in state. Although the dFINDE is designed to climinate temporal discretization errors in energy, it does not necessarily reduce those in state. In contrast, the Dormand- -ormand-Prince integrator is designed to suppress temporal discretization error in state.
21

---

Published as a conference paper at ICLR 2023
Therefore, there is no guarantee that the dFINDE improves the prediction performance when defined using errors in state. Conversely, the experimental results in Table 3 demonstrate that the dFINDE is superior to the base model and cFINDE in VPT. This is because dFINDE reduces the modeling errors rather. For the mass-spring system, the governing equation is already known as an ODE and is discretized by the dFINDE, leading temporal discretization errors. However, when dFINDE learns dynamics from data, the training data points are already sampled in discrete time, and the dFINDE predicts future states in discrete time. Therefore, no temporal discretization error occurs, and we obtain only the advantages of exactly preserving the first integral.
This type of paradox has been repeatedly discovered in previous studies. For example, the leapfrog integrator and discrete gradient method are second-order methods. However, they are superior to the Dormand-Prings integrator when combined with neural networks and learning dynamics from data (Matsubara et al., 2020). For better learning (i.e., smaller modeling errors), the preservation of specific properties of target systems more important than the order of accuracy.
D.2 SYMBOLIC REGRESSION OF FOUND FIRST INTEGRALS
Using gplearn (based on genetic programming), we performed symbolic regression of the first integrals V found by the neural network. We prepared addition, subtraction, multiplication, and division as candidate operations, used Pearson's correlation coefficient as the evaluation criterion, set the early stopping threshold to 0.9, and set the population size to 10,000, We set the other hyperparameters to their default values, c.g., the maximum number of generations was 20.
We summarize the regression results of the HNN with eFINDE for K = 2 trained using the two- body dataset in Table Al Note that Pearson's 's correlation coefficient is invariant to biases and scale factors. FINDE is also invariant because it only uses the directions of the gradients of first integrals. Hence, removed biases and scale factors from the regression results. When the focus is on the symbolic regression of the training data, V1, Vi, V22 and 1/2 for trials 0, 1, 2, and 3 are identical to the linear momentum in the x-direction up to scale factors; recall that we set mi m2 1.0 and see Eq. (A17). V2, V2, Vi, and V1 for trials O, I, 2, and 3 are also identical to the linear momentum in the g-direction. V and V/2 for trial 4 are weighted sums of the linear momenta in the x- and g-directions; in particular, they can be regarded as the linear momenta in the (1, -1)- and (1, directions, respectively.
When the quantities Vi (1)) and (u( are first integrals, any function of only V1 (u), V2(u),, and arbitrary constants is a first integral functionally dependent on V1 (u) and V2(u). Thus, in principle impossible to re-discover a first integral as well-known symbolic expression, and failure in symbolic regression is not a problem in any way. Previous studies introduced certain constraints (such as "gauge fixing") for symbolic regression (Liu & Tegmark 202T); combination of such method may improve the results. However, recent studies on neural networks have revealed that typical initialization and training procedures tend to learn simple functions (Barrett & Dherin [02T]] Cao et al. (022). Additionally, the symbolic regression limited the depth of the computation graph, biasing the results toward simple functions; hence, the found first integrals were identical to the well-known forms and were separated in the x- and y-directions in most cases.
The same is true for the symbolic regression of the test data, except for V1 for trial O, which had a small perturbation a. Because of the limited extrapolation ability, neural networks cannot always accurately represent functions outside the training data range. Once first integrals are found by FINDE and identified as equations by symbolic regression, one can use the equations instead of neural networks, ensuring the preservation of first integrals in the entire domain. From these results, we can conclude that cFINDE identified the linear momenta.
The state of the KdV dataset has 50 elements, which too large to apply a symbolic regression. For the 2-pend and FizZugh-Naguno datasets, we did not find consistent equations of first integrals. For example, the symbolic regression identified a quantity z / as a first integral in the 2-pend dataset, which is not directly related to well-known first integrals. When the angle 0, of the upper rod is small, y/1 takes a value close to -1, and the quantity of - y/1 is close to 2f + vf, which is a well- first integral, namely the square I3 of the upper rod length 41. It is difficult to determine whether this inaccuracy is because of the training of FINDE or symbolic regression. There may still be room for improvement in the training of FINDE symbolic regression.
22

---

Published as a conference paper at ICLR 2023
Table Al: Symbolic Regression of First Integrals Found in Two-Body Problem
Training Data
Test Data
Trial
1/1
V/2
0
Mai +/2
"y1 1/2
Mal 1u2 +
Vy1 + 1y2
Uyl 1x-
Uy1 +y2
Us +2²2
Vyl +uy
2
eyy I 14/2
Mal I Va2
Uyl I 1y2
Val I Uall
3
Pyl 1y2
Pz1 Px2
Myy +2y2
Fxrl +122
Myl + Uar Uy1 2y2
Vel + H/2 + Uyt 4y2
Mal + My2
Uyl
Vy2
Met + Ve2 + Vyy + +222
We removed biases and scale factors. a 0,003(y1
yi (M,2 1 1/ 32)
1.402).
Table A2: Results with Known Holonomic Constraints.
2-pend
2-body
Model
1-stepl
VPTV
1-stepl
VPT1
NODE
0.82 +0.0
0.110 20035
144.21 -1265
0.134 20.014
HNN Greydanus et al. 2019)
6220. 26 191.57
0.002 100000
5.17 -0.7
0.362 10.026
CHNN Finzi et al. 20206
0.07 10.00
0.928 20036
(not working)
FONDE
0.71 10.04
0.461 20071
163.64 79.79
0.147 20:024
HNN+FFINDE
236.51 #7.15
0.020 20002
8.32 -0.4
0.476 +0,400
D.3 COMPARISON WITH MODEL OF KNOWN HOLONOMIC CONSTRAINTS
The double pendulum (2-pend) is classified as a constrained Hamiltonian system. CHNN was pro- posed for cases when holonomic constraints known (Finzi al. 2020b). We evaluated compari- son methods under the assumption that the holonomic constraints were known. We summarized the results in Table A2 The HNN, without constraints, completely failed learn the dynamics. This is unsurprising because the dynamics of the double pendulum is outside the scope of the HNN. The two known holonomic constraints lead to two constraints involving the velocity; the CHNN took into account all four known constraints and worked remarkably. The HNN with cFINDE was given all four known constraints as the first integrals, but did not work properly. The original purpose of projection methods is to eliminate temporal discretization errors of first integrals but not to change the class to which the dynamics belong. Therefore, when a target system not a subject of the base model. the base model with FINDE does not work. The NODE learns an ODE in general way, and thus constrained Hamiltonian systems are included in its subjects. Given all four known constraints, the NODE with cFINDE worked better but never surpassed the CHNN.
However, the CHNN works only for Hamiltonian systems in the canonical form with holonomic constraints. We also evaluated comparison methods using the 2-body dataset under the assumption that the linear momenta were known as first integrals. The CHNN attempted to obtain the inverse of singular matrix and could not learn the dynamics. In contrast, the cFINDE improved the perfor- mances of both NODE and HNN.
Existing methods (e.g., HNN and CHNN) assume geometric structures (e.g., Hamiltonian struc- ture) described in Appendix A in order to guarantee conservation laws. When multiple structures are assumed at the same time, they must be integrated using appropriate prior knowledge. If it is possible, it would achieve extremely high performance. Otherwise, the geometric structures would conffict with each other and would not produce an appropriate model. This the reason why CHNN failed to learn the 2-body dataset and HNN+FINDE failed to learn the 2-pend dataset. In contrast, NODE+FINDE does not assume any geometric structure and assumes first integrals in the most general way, being available to any situation. Hence, FINDE can assume or more first integrals without changing anything.
When the detailed properties of target systems are known, one can choose the best models. If the chosen model is inappropriate, the training procedure totally fails. FINDE provides a better alternative when prior knowledge is limited. Moreover, constrained Hamiltonian system can have first integrals other than holonomic constraints and the Hamiltonian. In this case, the CHNN with FINDE is potentially the best choice.
23

---

Published as a conference paper at ICLR 2023
Table A3: Results of NODE with cFINDE on Training Set of 2-Pend Dataset.
2-pend
Model
K
1-stepl
VPTT
NODE
0.76 =0:02
0,966 =0:007
1 0.72
-0066
0,974
2
0.69
0,981 =0.014
cFINDE
3
0.63
0,994 00.002
4
0.67
20055
0.990 20.005
5
0,998 +0.000
6 9.93
0,126 10,00
D.4 REASON FOR HIGH PERFORMANCE AND How TO DETERMINE NUMBER OF FIRST INTEGRALS
The theoretical explanation for the high performance of neural networks (e.g.. HNN) that assume first integrals for physical phenomena an open question. [annaii et al. (202T) has theoretically shown that neural networks (c.g., CNNs and GNNs) with symmetry have faster learning conver- gence, and we consider this approach can be applied to the above question. At least for cFINDE and dFinde, we have an intuitive but not rigorous explanation; assuming one more first integral (i.e., increasing K by reduces the number of degrees of freedom in the dynamics by 1, narrows the hypothesis space, accelerates learning convergence, and suppresses generalization errors.
As shown in Table the performance of cFINDE and dFINDE sensitive to the assumed number K of first integrals. Because K is a hyperparameter, it is basically a subject to be adjusted through evaluations on a validation set. With inappropriately large K. both cFINDE and dFINDE dropped their performance significantly. See the results of the 2-pend and FitzHugh-Naguna datasets for K = and K = 3, respectively.
However, the performance drop be found even with the training set. Table A3 summarizes the prediction performance on the training set of the 2-pend dataset. As was the case with the test set, the performance significantly dropped at K 6. This is because NODE with cFINDE for K assumes the submanifold M' to be 2-dimensional. The submanifold M' is in fact 3-dimensional, so NODE with cFINDE for K = 6 is incapable of learning the dynamics and performs poorly even on the training set. Hence, the training set is enough to avoid a fatally inappropriate K,
Alternatively, K can be determined by using other methods (c.g., .ukunaaa & Olsen (197] Liu & Tegmark (202T)). Although these methods have some drawbacks introduced in Appendix they may be complementary FINDE.
D.5 COMPARISON WITH MODIFIED NEURAL PROJECTION METHOD
The neural projection method (NPM) also employs projection method Yang et al. 2020). Using a manner similar Newton's 's method, it enforces the constraint C'uu) by the projection of the state u under the assumption that the quantity C(u) is always zero. This assumption holds for some cases (e.g., holonomic constraints in a fixed environment), but not for most first integrals, whose values depend on initial conditions.
For example, the linear momentum in the z-direction of the two-body problem is the first inte- gral expressed V(u) = myvzi (() + (l).,((). This quantity V is constant within trial (i.e., V(u(t)) = (O)) and varies between trials depending on the initial speed Uy1 (0) and (0) The total energy, the total mass, and many other first integrals depend on the initial condition in the same manner; hence, they are outside the scope of the NPM. In contrast, by imposing the constraint on the gradient VV 0 or discrete gradient VV 0, our proposed FINDE keeps the quantity V constant and can handle any first integrals.
For comparison, replaced the constraint C(u) 0 with C(u+11 w²) V( V(u")= and adopted the NPM to first integrals varying from trial to trial. We evaluated the modified NPM using the 2-pend dataset. Because the modified NPM is discrete-time projection method,
24

---

Published as a conference paper at ICLR 2023
Table A4: Comparison with Neural Projection Method (NPM)
K
dFINDE (proposed)
modified NPM
1-step/
VPT
1-stepl
VPT1
successful
0.75 20.10
0.152 20,07
0.73 20.8
0,150 20014
5/5
2
0.74 20.05
0.271 H0.111
0/5
3
0.69 10.05
0.447 10.081
(0.69 10.00)
(0.138 20.00)
1/5
0.71 2003
0.454 +0,000
(0.72 2000)
(0.383 20023)
0.86 40.09
0,591 +0.057
40.11
0,364 #0.134
58.88 =22.98
0.037 10,039
(1.29 10.20)
(0.103 20016)
3/5
we compared it with discrete-time version of the proposed FINDE (dFINDE). The results are summarized in Table A4
The dFINDE successfully learned the dynamics in all trials, but the modified NPM failed to learn the dynamics in half the trials (see the rightmost column for the numbers of successful trials out of 5). The modified NPM often encountered of the underflow of the time-step size or a division by the zero gradient of the first integral. Even when the learning was inferior to that of the dFINDE. The modified NPM solved the optimization problem in Eq. 3 at every step, but it sometimes diverged or failed converge, especially the early phase of learning. The NPM was successful for fixed environments but might be unsuited for general first integrals varying from trial trial. However, the dFINDE does not require solving an optimization problem during training, making the learning process robust against randomness such as initialization.
successful, the performance of the NPM
25


Reviews:
[{"review_id": "cNOMq8X16r6", "reviewer": "Reviewer_vr6Y", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper introduces a new method to preserver first integrals in neural network modeling called first integral-preserving neural differential equation (FINDE) with the claim that FINDE can preserve and discover first integral in data without needing explicit information about the underlying structure. The authors propose two versions of FINDE: cFINDE (continuous integration) and dFINDE (discrete integration). In its introduction, the paper outlines the importance and challenges related to finding first integrals from data and which first integrals can be especially relevant for computer simulations (Energy, momentum, mass, etc). Next, the authors formally define first integrals and describe related work of first integrals in numerical analysis and neural networks that preserve first integrals, including Neural Ordinary Differential Equations, as well as works that have built on top of it and different kind of architectures (holonomic neural networks - HNN). The paper then introduces and defines the two versions of FINDE: cFINDE and dFINDE with relevant theorems explaining important preservation properties related to first integrals for both formulations.\n\nIn their experiments, the authors evaluate cFINDE and dFINDE on relevant datasets that contain different kinds of first integrals. In their implementation, the authors leverage an HNN,  apply related loss formulations to train cFINDE and dFINDE, and define valid prediction time (VPT) as an evaluation metric. The first set of experiments investigate whether first integrals are preserved in FINDE via a mass spring example. The second set of experiments show results on finding known integrals with FINDE compared to using vanilla HNN with the results indicating that FINDE methods generally perform better. The third set of experiments deal with finding first integrals in unknown systems using a set of examples (KdV, pendulum). The results generally indicate that FINDE methods are able to reproduce the ground truth compared to neural ODEs.", "summary_of_the_review": "The paper introduces a novel method and provides a relevant set of experiments to showcase the abilities of the method. Improvement can be made to further strengthen the paper in subsequent revisions.", "strength_and_weaknesses": "**Strengths**\n\n* The paper introduces a novel formulation for a relevant problem in applying neural network to physical systems and computer simulations thereof. The FINDE methods are described and derived with detail and their properties tested as part of the experiments.\n* The paper provides a relevant set of experiments to assess the properties and performance of FINDE, both cFINDE and dFINDE and seems to achieves superior performance compared to other methods.\n\n**Weaknesses**\n\n* The paper could be improved by having more detailed Figure captions across all figures, especially the one outlining the results (e.g. Figure 1, Figure 4 and Figure 5). More detailed captions would provide greater clarity on what the data show and how to interpret the results.\n* The authors could further clarify potential sources of error in FINDE. The current version explains how dFinde resolve the time discretization error in cFINDE, but it would be good to have a more information on other sources of error created by the approximation and they can potentially be resolved. E.g. In any finite difference scheme there is approximation error, such as truncation errors coming from the discretization of the equation.\n\n**Additional Questions**\n\n* You mention that FINDE preserves multiple quantities at the end of Section 2. Could you explain more how that occurs? Specifically how does preserving multiple quantities (e.g. energy and mass) differ from preserving just one quantity (energy) related to your method?\n* Given the additional advantages of dFINDE, in what situations would someone choose cFINDE over dFINDE?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "MbCqgjZtxA", "reviewer": "Reviewer_pYdf", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper asks a very pertinent question as to how we can incorporate knowledge of conserved quantities when solving PDEs via nets. The suggested method is quite unclear. ", "summary_of_the_review": "Given the weaknesses described above, I feel there are major gaps here - either very fundamental issues have gotten overlooked from the study or the presentation has failed to communicate them to the reader.  The paper needs a major rewrite in the very least.  ", "strength_and_weaknesses": "The main strength of the paper is that it raises and important question and that it has tried some serious benchmark tests. Also interestingly this paper tries to marry the ideas of symbolic regression with ideas of neural training to discover conserved quantities. Definitely a worthwhile attempt. \n\nPrima facie there seem to be major weaknesses in the paper.  \n\nFor a start V seems to be a neural net whose output coordinates encode the conserved quantities. But is V being trained? If yes, how? If I understand correctly -- the function f in Theorem 1 is also a neural net. Is that being trained? Its entirely unclear to me as to how these nets are getting determined. If there is a loss function somewhere that is getting minimized then that is extremely unclear. The paper is in dire need of a pseudocode! \n\nThe testing of preservation of the conserved quantities looks very ambiguous to me. While the authors want to use this as one of the main tests it doesn't seem to have a clear definition till Appendix A (Equation A13?) - and even then it seems to be a function of time while the Table 3 seems to suggest that this is a number that characterizes the whole solution.  And if we focus on Table 3 then it seems that at no value of K is the new method beating the benchmarks in this 1-step criteria for this 2-body example. Isn't that a strong argument against the new method? \n\nFurther, in Appendix D.2 the authors seem to clearly state that their approach of symbolic regression is not leading to a correct identification of the conserved quantities for their 2-pendulum and the FitzHugh-Nagumo datasets. So what exactly is the claim of the paper at all? Doesn't this beat the central claim in the paper? \n\nUsing KdV as a test-case is a quite tricky. Firstly there seems to have be a typo in equation A18 when quoting the KdV. Further KdV has infinite number of conserved quantities. The authors claim to have cut that down to some finite number by some discretization. This is entirely unexplained in the paper - and potentially a very fatal issue since that sounds like trying to entirely change the target from trying to solve KdV to something else! ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "8iPqdf-Fth", "reviewer": "Reviewer_RgJp", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "This paper proposes first integral-preserving neural differential equation (FINDE) that can learn a dynamical system and its invariant quantities (i.e., first integrals) from observed data *without* prior knowledge (e.g., conservations of energy and momentum) or assumed geometric structures (e.g., Hamiltonian or Lagrangian) on the target system. FINDE consists of two neural networks; a neural network $\\hat{f}(u) \\in \\mathbb{R}^N$ for modeling the target dynamics (i.e., neural ODE) and another neural network $V(u) \\in \\mathbb{R}^K$ that outputs $K$ first integrals. Then, FINDE projects the neural vector field $\\hat{f}(u)$ to the directions in which $V(u)$ is preserved, by using the projection method. Thus, it guarantees the state $u$ evolves on $N-K$ dimensional sub-manifold (level set) which preserves first integrals. The optimization procedure of the projection method is well-formulated by using Lagrange multipliers, and it yields a novel type of neural ODE (continuous FINDE, cFINDE) that can learn conservative dynamics with repsect to the modeled frist integrals. The authors also present a discrete version of the proposed FINDE (dFINDE) to preserve the learned first integral under the discretization of numerical integrators. The authors theoretically show cFINDE and dFINDE exactly preserve the modeled first integrals $V(u)$ under the continuous and dicrete cases, repectively. The authors validate their approach with four dynamical systems including simple two-body problem and complicated KdV system.", "summary_of_the_review": "Overall, I believe the paper is very interesting, tackles an important problem, and thus is worthy of publication. However, the concern on $K$ tuning prevents me from recommending a clear accpetance of the paper in its current form. Currently, my evaluation is weak accept, however, I would like to increase my score if my concern is well addressed.", "strength_and_weaknesses": "Strength:\n\n1. The motivation of this paper is clear. Extracting underlying invariant quantities as well as learning dynamics *without* prior knoweldge on the target system is definitely important and challenging problem in physics. This paper provides a solid contribution to such an important problem. It might be one of good references for automated scientific discovery via machine learning tools.\n\n2. The paper is technically sound. The used projection method seems to be a principled way to preserve arbitrary first integrals. The authors properly modifies the projection method optimization as a learnable model (FINDE) based on neural ODE framework. \n\n3. The paper is clearly written and easy-to-follow.\n\n4. This work is highly reproducible. All used model architecture and hyper-parameters are well summarized in the main text of this paper. Experimental codes were made available.\n\nWeakness:\n1.  $K$, the output dimension of $V(u)$, is very important hyper-parameter because it directly determines the number of invariant quantities. Inappropriate setting of $K$ rather deteriorates the predictive performance of the model as depicted in the experimental results. In the current version of this paper, the authors do not clearly mention how one can determine such a critical hyper-parameter for unknown dynamical systems. It will be nice if the authors elarborate it.\n\n2. I found a prior work [1] that seems to be very relevant to the proposed method. [1] uses PCA to find $N-K$ dimensional sub-manifold of trajectories and extracts $K$ invariant quantities. It will be nice if the authors cite [1] and compare pros and cons of the proposed method with [1] briefly.\n\n***\n[1] Liu, Z., & Tegmark, M. (2021). Machine learning conservation laws from trajectories. Physical Review Letters, 126(18), 180604.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "5PgnAW9Rjq", "reviewer": "Reviewer_USE3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "This paper proposes a neural network based approach to preserve and discover first integrals in the underlying target systems. In this paper, authors propose two instances of the approach, namely cFINDE and dFINDE, which works in continuous-time and discrete-time respectively. The proposed FINDE method is a great exploration in the area of preserving and/or discovering invariant quantities when learning dynamical systems with deep learning methods, and it is a generalization of existing deep learning based methods. As I see, the FINDE method makes a solid contribution.  ", "summary_of_the_review": "As stated in the above comment sections, this paper has a solid contribution with a novel approach. It is a good paper.", "strength_and_weaknesses": "The main strengths of this paper are \n\n$\\mathbf{1}.$ on the proposal of enhancing the neural network based approach to learn dynamical systems by proposing an additional structure/equation to preserve and/or discover invariant quantities. \n\n$\\mathbf{2}.$ This paper is very well organized, written and illustrated with graphic results. The background knowledge is adequate for readers (relevant in deep learning fields) to understand the rationale and the intuition driving the proposal, and the paper is very readable.  \n\nThere are a few weakness, which could make the paper better but should not overshadow the strength\n\n$\\mathbf{1}.$ The numerical results can be explained and presented better. For example, the switch between dFINDE and cFINDE could be explained better as they often make a bit confusion in the text; the plots can be improved with better presentation: for instance, Figure 3 - it can be shown as four subplots, each of which can compare HNN vs HNN+cFINDE; the plots could confuse audience a bit: for instance, it is hard to see the ground truth. \n\n$\\mathbf{2}.$ It is very interesting to see the big difference induced by using different values of $K$. As discovering the invariant quantities is one of the cores of this approach and also as author addressed in the conclusion, conducting more comprehensive and in-depth numeric investigation on the role of $K$, with more detailed discussion, along with more analysis in theory, would benefit a lot in helping audience understand and the field leverage the proposed methods.    ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster