KAM Theory Meets Statistical Learning Theory: Hamiltonian Neural Networks with Non-Zero Training Loss

Many physical phenomena are described by Hamiltonian mechanics using an energy function (the Hamiltonian). Recently, the Hamiltonian neural network, which approximates the Hamiltonian as a neural network, and its extensions have attracted much attention. This is a very powerful method, but its use in theoretical studies remains limited. In this study, by combining the statistical learning theory and Kolmogorov-Arnold-Moser (KAM) theory, we provide a theoretical analysis of the behavior of Hamiltonian neural networks when the learning error is not completely zero. A Hamiltonian neural network with non-zero errors can be considered as a perturbation from the true dynamics, and the perturbation theory of the Hamilton equation is widely known as the KAM theory. To apply the KAM theory, we provide a generalization error bound for Hamiltonian neural networks by deriving an estimate of the covering number of the gradient of the multi-layer perceptron, which is the key ingredient of the model. This error bound gives an $L^\infty$ bound on the Hamiltonian that is required in the application of the KAM theory.


Introduction
Many physical phenomena are described by energy-based theories, such as Hamiltonian mechanics (e.g., Furihata and Matsuo (2010) where u : t ∈ R → u(t) ∈ R N , H : u ∈ R N → H(u) ∈ R, and S is a skew-symmetric matrix. H represents the energy function of the system. In recent years, there has been a lot of research on predicting the corresponding physical phenomena by learning the energy function H in such equations with a neural network H NN (e.g., Chen et al. (2020); Cranmer et al. (2020); Greydanus, Dzamba, and Yosinski (2019); Matsubara, Ishikawa, and Yaguchi (2020); Zhong, Dey, and Chakraborty (2020)); however, to the best of our knowledge, theoretical analysis of such models has not been performed sufficiently, except for SympNet (Jin et al. 2020b) for the Hamilton equation, where the universal approximation theorems for discrete-time neural network models are provided.
In this paper, we focus on theories of the properties of the most fundamental model, comprising Hamiltonian neural networks (HNNs) (Greydanus, Dzamba, and Yosinski 2019) and their extensions in practical situations, where the learning error is not completely zero. In this case, the trained model can be regarded as a perturbed Hamiltonian system due to the modeling error of the energy function. In addition, S is a general skew-symmetric matrix and hence (2) can model Hamiltonian partial differential equations (Matsubara, Ishikawa, and Yaguchi 2020). In mathematical physics, perturbed Hamiltonian systems are well studied. For example, because whether the solar system will continue to exist in the future is of fundamental interest in astronomy, the stability of the solar system under perturbations is a very important issue that has been studied for a long time (e.g., Féjoz (2013); Laskar (1996)). The Kolmogorov-Arnold-Moser (KAM) theory gives an answer to questions of this type; essentially, periodic motions of such systems are stable under small perturbations. The stability of periodic motions is of particular importance in science. In addition to the stability of celestial systems, the recursive nature of physical phenomena is also of interest in physics. For example, in the famous numerical experiments using the Korteweg-De Vries (KdV) equation by Zabusky and Kruskal (1965), it was confirmed that the waveform given as the initial condition initially collapsed due to severe oscillations, and then returned to its original shape after a long time. Whether such phenomena can be reproduced by deep learning models is an important problem that greatly affects the usefulness of deep physical models.
In this paper, we give an answer to this question by combining the KAM theory and statistical learning theory. Because trained models are also perturbed Hamiltonian systems, it is expected that the periodic behaviors of the systems are preserved even if the loss function does not vanish completely provided it is sufficiently small. However, this expectation cannot be proved in a straightforward way, because the application of the KAM theory requires that the energy function be close enough to that of the true dynamics in the whole phase space. Noticing that the error of the Figure 1: Outline of our main theorem. The first main result is the generalization error of the energy function, which is proved by the statistical learning theory. By combining the first result and the KAM theorem, we prove that the quasi-periodic behaviors of target systems are preserved. energy function in the whole space is essentially the generalization error, we overcome this difficulty by combining the KAM theory with statistical learning theory. This illustrates that combinations of statistical learning theory and dynamical systems theory can lead to powerful results. Indeed, combinations of this kind may be applicable to other stability results in dynamical system contexts, for example, the stability of solitary wave solutions.
Importantly, for the neural network models to be close to the true dynamics, we need a universal approximation theorem and also a generalization error bound. In this paper, we also provide such results for HNNs.
Regarding the generalization error bound, because the derivative of a multi-layer perceptron is used in HNNs, a bound for the derivative is required. To this end, we estimated the covering number of the derivative of multi-layer perceptrons. An L ∞ bound on the error in the Hamiltonian is also provided, which is required for application of the KAM theory.
In addition, we show a universal approximation theorem for a model with the coordinate transformation This model is indispensable in practice; to apply HNNs, the data are given in the canonical coordinate system because the equation of motion is in the form of the Hamilton equation (1) only when the state variables are represented by the canonical coordinate system. However, this coordinate system depends on an unknown Lagrangian and hence the energy function. Hence, the coordinate system must be also learned from data by using, for example, neural networks.
In addition, this model also represents other energy-based physical models beyond the Hamilton equation. See the Appendix for details. The main contributions of this paper are as follows.
1. Combination of the KAM theory and statistical learning theory for HNNs with non-zero training loss to prove the existence of quasi-periodic behaviors (see Figure 1). 2. Derivation of a generalization error bound for HNNs.

Development of a universal approximation theorem
for HNNs and other energy-based physical models with coordinate transformations.

Related work
Many studies of neural network models for phenomena that can be modeled by energy-based equation (1) have been put forward. Among them, the most basic studies are neural ordinary differential equations (Chen et al. 2018) and HNNs (Greydanus, Dzamba, and Yosinski 2019). In particular, extensions of HNNs have been intensively developed. Describing them all is beyond the scope of this paper, but some examples are given here. In Toth et al. (2019) HNNs were extended to latent variable models. Other studies, such as DiPietro, Xiong, and Zhu (2020); Xiong et al. (2021);Zhong, Dey, and Chakraborty (2020), focused on the symplectic structure of the Hamilton equation. For Noether's theorem, which is a fundamental theorem in classical mechanics, several studies (Bharadwaj, Li, and Demanet 2020;Bondesan and Lamacraft 2019;Finzi et al. 2020) developed methods related to symmetry and conservation laws. In addition, a discrete-time model that preserves the energy behaviors was constructed in Matsubara, Ishikawa, and Yaguchi (2020). In Galioto and Gorodetsky (2020), HNNs were combined with a Bayesian approach.
Methods applied to the framework of classical mechanics other than Hamiltonian mechanics include those in Cranmer et al. (2020); Desai and Roberts (2020); Saemundsson et al. (2019), which are methods for Lagrangian formalism. In Jin, Lin, and Li (2020), reinforcement learning was applied to the variational principle. A simplified model formed by introducing constraints was proposed in Finzi, Wang, and Wilson (2020). In Jin et al. (2020a), HNNs were extended to the Poisson system, which is a wider class of mechanical equations. There are also a number of proposals that integrate them with more advanced deep learning techniques, namely, graph networks (Sanchez-Gonzalez et al. 2019), recurrent neural networks (Chen et al. 2020), and normalizing flows . As an application-oriented approach, Feng et al. (2020) designed a microscopic model for structural analysis.
However, as far as the authors know, there is no theoretical research other than the universal approximation theorems for Hamiltonian mechanics in SympNet (Jin et al. 2020b), in which a certain kind of neural network is shown to have universal approximation properties for symplectic maps. The difference between their results and ours is that (1) we analyze the behaviors of a HNN with non-zero training loss by a combination of the KAM theory and statistical learning theory, (2) we provide a generalization error bound for NHHs, and (3) the universal approximation theorems in Jin et al. (2020b) are for discrete-time models, while ours are for continuous-time models.
Meanwhile, as an existing energy-based model, Hopfield neural network is known. Both Hopfield neural network and Hamiltonian network are derived from energy-based theories, and their dynamics are described by (1). Hamiltonian neural network is associated with a skew-symmetric matrix S and is a model of an energy-preserving, continuous-time, and deterministic physics phenomenon. Its output is the time-series of the state. Hopfield network is associated with a negative definite matrix S and exhibits a dynamics which is often energy-dissipating, discrete-time, and stochastic. It is a machine learning tool rather than a physical model, and its equilibrium point is treated as its output. Because their outputs are different, their theoretical properties should be discussed separately.

Brief introduction to Hamiltonian systems and the KAM theory
We briefly introduce some properties of Hamiltonian systems. Theorem 1 (Darboux). By an appropriate coordinate transformation, the matrix S is transformed into the normal form Definition 1. The function ω : is called the symplectic form. Using the symplectic form associates a vector field X F with each function F : R N → R by requiring For two functions F, G, the following operation is called the Poisson bracket: (4) Definition 2. A Hamiltonian system for which the state variable is N = 2M dimensional is integrable in the sense of Liouville if this Hamiltonian system has the first integrals (i.e., conserved quantities) F 1 , F 2 , . . . , F M with ∇F 1 (u), ∇F 2 (u), . . . , ∇F M (u) independent at each u and for all i, j: {F i , F j } = 0. For integrable systems, Theorem 2 is known. Theorem 2 (Liouville-Arnold). Suppose that for an integrable Hamiltonian system, constants c 1 , . . . , c M exist such is connected and compact. Then, there exists a neighborhood N comprising K, U ⊂ R n and a coordinate transform φ : (θ, J) ∈ T n × U → φ(θ, J) ∈ N (5) such that the transformed system is the Hamilton equation of which Hamiltonian H • φ depends only on J.
The variables J and θ are called action-angle variables. Theorems 1 and 2 roughly mean that integrable Hamiltonian systems can be written in the following form: Further, becauseH = H • φ depends on I only, it holds that d dt This shows that J is constant, and hence θ moves on the torus at a constant velocity. Because the velocities are typically not co-related to each other, the dynamics are "quasiperiodic." See, for example, Scott Dumas (2014) for more details.
As seen above, integrable Hamiltonian systems are quasiperiodic. Note that general Hamiltonian systems are not necessarily quasi-periodic and neither are HNNs. However, for the HNNs that are trained to model integrable systems, the quasi-periodic behaviors are preferably maintained. When the modeling error is sufficiently small, this is considered as a perturbation problem. The perturbation theory of Hamiltonian systems has been investigated from various perspectives. For example, perturbed integrable Hamiltonian systems are in general no longer integrable; hence, approximation of integrable Hamiltonian systems by integrable neural network models appears to be difficult. Fortunately, however, the KAM theory shows that even though the perturbed system is not integrable, it maintains the quasi-periodic behaviors described above under certain conditions.
The KAM theorem has many variants under various conditions. The following variant (Scott Dumas 2014) is typical: Theorem 3 (KAM Theorem). Let θ and J be the actionangle variables for a C ∞ integrable Hamiltonian H 0 : for the perturbed system H(θ, J) = H 0 (J) + εF (θ, J, ε) by F ∈ C ∞ , there exists ε 0 such that if εF < ε 0 , there exists a set of M -dimensional tori that are invariant under the perturbed flow. For each invariant torus, the flow of the perturbed system H is quasi-periodic. In addition, the set of invariant tori is large in the sense that its measure becomes full as ε → 0.
Remark 1. The last sentence -the set of invariant tori is large in the sense that its measure becomes full as ε → 0corresponds to the non-existence of so-called resonance. If the perturbation added to the system is in resonance with the original system, the perturbation may grow rapidly and the behavior of the system may change significantly. This statement assures that for small perturbations, such a situation almost never occurs.
Remark 2. It may be difficult to check whether the target system is integrable by using given data. One possibility is application of the Koopman operator, which makes it possible to find the conserved quantities that the given data may admit. If a sufficient number of conserved quantities exist, it is highly likely that the target system is integrable.

Universal Approximation Properties of HNNs
For HNNs to be close to the true dynamics, a universal approximation theorem and a generalization error analysis are needed. First, we show universal approximation theorems. We first define some notation to describe the theorem. C m (X) with the topology of the Sobolev space W p,m (X) is denoted by S m p (X), where W p,m (X) is a space of functions that admit weak derivatives up to the mth order that bounds L p -norms. Hence, S m p (X) is the space of functions in W p,m (X) with (usual) derivatives; for details on the Sobolev theory, see Adams and Fournier (2003). L p -norms of functions are denoted by · L p , and those of vectors by · p .

Universal approximation theorem for HNNs
The following theorem shows the universal approximation property of HNNs.
Theorem 4. Let H : R N → R be an energy function with the equation where u : t ∈ R → u(t) ∈ R N and S is a skew-symmetric N × N matrix. Suppose that the phase space K of this system is compact and the right-hand side S∂H/∂u is Lipschitz continuous. If the activation function σ = 0 belongs to S 1 2 (R), then for any ε > 0, there exists a neural network H NN for which holds. In addition, if the energy function is C ∞ , the function can be approximated by a C ∞ neural network provided that the activation function is sufficiently smooth. An outline of the proof is as follows. From the existence theorem of a solution to ordinary differential equations, the right-hand side of the equation must be Lipschitz continuous to guarantee the existence of a solution. Therefore, the differential of the H is Lipschitz continuous, which implies the smoothness of H. Then, the approximation property is obtained from the universal approximation theorem for smooth functions by Hornik, Stinchcombe, and White (1990). For the detailed proof, see the Appendix.
HNNs with a coordinate transformation The practical use of HNNs is hampered by the fact that the state variables must be represented by a specific coordinate, such as the generalized momentum; however, the derivation of the generalized momentum requires the energy function, which is unknown. For example, the double pendulum in Figure  2 exhibits the dynamics shown in Figure 3. These are predicted by the models that are trained from the data of the state variables and their derivatives, not those of the generalized momenta. HNNs failed to solve such problems because the data were not given in the canonical coordinate system. See the Appendix for details. Based on this, we here propose a model with a coordinate transformation, such as the transformations that appear, for example, in Rana et al. (2020); Jin et al. (2020a).
Suppose that, although the given data point x(t) is not represented by the canonical coordinate system, the data point x(t) can be transformed into the canonical coordinate system by an unknown transformation We show that model (7) admits the same energetic property as the original equation and also the universal approximation property. Theorem 5. The model (7) admits the energy conservation law in the sense that dH NN /dt = 0.
Theorem 6. Let H : R N → R be an energy function for the equation Suppose that the phase space K of this system is compact and the right-hand side ∂H/∂u is Lipschitz continuous. Suppose also that u is a C 1diffeomorphism. If the functions σ = 0 and ρ = 0 belong to S 1 2 (R), then for any ε > 0, there exist neural networks H NN with the activation functions σ and u NN with ρ for which For proofs, see Appendix.

Generalization Error Analysis of HNNs
Next, we derive a generalization error bound for the standard HNN (2) by employing a technique from statistical learning theory. More precisely, we adjust the technique so that an estimate on the energy gradient can be obtained. Remark 3. Although the standard HNN without the coordinate transformations is considered below, the results can be extended to the general energy-based model with the coordinate transformations if the matrix (∂u/∂x) −1 is bounded.
Lemma 1. Let X and Y be arbitrary spaces, F ⊂ {f : X → Y } be a hypotheses class, and L : Then, for any δ > 0 and any probability measure P , we obtain, with a probability of at least (1 − δ) with respect to the repeated sampling of P n -distributed training data, the following: The Rademacher complexity is known to be bounded by using the covering number.
Thus, if the covering number is estimated for a HNN, the bound on the generalization error is obtained. To this end, we suppose that the model is trained by minimizing the pnorm of the error in the right-hand side of the model. More precisely, for the hypothesis h : u j → S ∂HNN(uj ) ∂u , we consider the loss function where u i are training data. We denote the Lipschitz constant of the loss function associated with the above by ρ p . Of course, p = 2 is typically used; however, we show below that p > 2M is useful to obtain an L ∞ bound on the Hamiltonian.
Remark 4. The training can be performed also by using the symplectic gradient: In that case, the results will be slightly modified using the norm of S; however, we omit this for simplicity.
A bound of the covering number is derived as follows.
Theorem 7. Suppose that the hypotheses class F consists of multi-layer perceptrons f NN that have ρ σj -Lispchitz activation functions σ j (j = 1, . . . , n l ), for which the derivatives are ρ j -Lipschitz continuous and bounded by sup |σ j | < c σj . Suppose also that the matrices A j (j = 1, . . . , n l + 1) in the linear layers in the perceptrons have the bounded norm |A j | < c Aj : with the ρ p -Lipschitz continuous loss function L. In addition, suppose that the phase space is compact so that the data u i (i = 1, . . . , n) are in a bounded set with the bound u i < c u . Then, the covering number of G is estimated by To prove this theorem, we use the following lemmas, which are typically used to estimate the covering numbers (Shalev-Shwartz and Ben-David 2014).
Lemma 4. Suppose that functions φ i : R → R, i = 1, 2, . . . , n are ρ-Lipschitz continuous. Then, for V ⊂ V n , Proof of Theorem 7. To simplify the discussion, we will estimate the covering number of the following perceptron Because the proof for general cases is exactly the same, we need to estimate the covering number of the gradient of f NN , which is written as where Dσ 2 and Dσ 1 are Jacobian matrices. These Jacobian matrices are evaluated at u = A 2 σ 1 (A 1 u + b 1 ) + b 2 and A 1 u + b 1 , respectively. We first estimate the covering number associated with Dσ 2 . Dσ 2 has the same architecture as a multi-layer perceptron, except that the last activation function is replaced by the differential σ 2 of σ 2 : Assuming that the input data are in the ball B cu with radius c u , we obtain Then, because the norms of A 1 , A 2 are bounded by c A1 , c A2 , matrix multiplications by these matrices are c A1 -and c A2 -Lipschitz continuous, respectively. In addition, σ 1 is ρ 1 -Lipschitz and σ 2 is ρ 2 -Lipschitz continuous. Therefore, the covering number associated with Dσ 2 is estimated by Finally, because σ 1 is assumed to be bounded by c σ1 , the norms of the matrices other than Dσ 2 in ∇f NN are bounded as follows: Because the loss function is assumed to be ρ p -Lipschitz continuous, we obtain the estimate L ∞ estimate on the error in the Hamiltonian The generalization error analysis in Theorem 7 shows that, at a certain probability, the expectation of the loss function can be bounded. If this bound certainty holds and if the training is performed by minimizing the p-norm with p > 2M , we can derive an L ∞ estimate on the Hamiltonian for the standard HNN (2) by applying the Poincaré inequality and the Sobolev inequality under Assumption 1.
Assumption 1 There exists a density f P for measure P with inf f P > 0. Remark 5. The condition p > 2M is not required in practice because of the well-known equivalence of the norms in finite dimensional spaces; for example, if the standard 2norm is small enough, then the p-norm is also small. However, when the dimension 2M is large, the 2-norm needs to be very small to bound the p-norm because the constant in the inequality used to bound the p-norm depends on the dimension. Therefore, it is preferable to minimize the p-norm in such cases. Theorem 8 (Poincaré inequality). Suppose that 1 ≤ p ≤ ∞ and Ω ⊂ R 2M is bounded. Then there exists a constant c p such that, for any H ∈ S 1 p (Ω), The constant c p is called the Poincaré constant. Theorem 9 (Sobolev inequalities). There exist constants c 1 , c 2 such that, if lp > 2M , By using these inequalities along with the invariance of the Hamilton equation under the constant shift of the energy function, we obtain an error bound on the Hamiltonian. From the above estimate, we obtain By using the density f P for the measure P , we obtain which gives us We note that the left-hand side is the Sobolev norm of the error in W p,l ; then, under the assumption that p > 2M , we can use the Sobolev inequality to obtain which ensure that H and H NN are close in terms of the function values.

KAM Theory for HNNs
The universal approximation property shown in the previous sections guarantees that the value of MSE can be made arbitrarily small by training; however, in actual training, a finite error remains. In this section, as an application of the generalization bound, we apply the KAM theory to theoretically investigate the trained standard HNN model (2) in such cases by assuming that the target system is integrable.
We make a few assumptions that are needed for the application of the KAM theory. Assumption 2 The dimension of the phase space is assumed to be 2M with M ≥ 2. Assumption 3 The target system is an integrable Hamiltonian system with the conserved quantities F 1 , . . . , F M . The series c 1 , . . . , c M exists such that K = ∩ M i=1 F −1 i (c i ) is connected and compact.
Under the above assumptions, from the Liouville-Arnold theorem there exist a neighborhood N of K, U ⊂ R n and a coordinate transform such that the transformed system is the Hamilton equation.
Following the usual setting of the KAM theorem, we consider the target system and the Hamiltonian equation in the transformed coordinate T n × U. Assumption 4 The Hamiltonian H : T n ×U → R of the target system is C ∞ and non-degenerate. The activation functions of the HNNs used are in C ∞ . Assumption 5 From the generalization error analysis in the previous section, we have essentially shown that if p > 2M , with at least probability 1 − δ, it holds that sup |H(u) − H NN (u)| < c 1 L train + c 2 R n + c 3 ln 1 δ n with constants c 1 , c 2 , and c 3 , where R n is a bound on the Rademacher complexity. We assume that the training was performed with p > 2M and the above statement certainly holds.
Using these assumptions, we obtain Theorem 10. Theorem 10. Let the threshold of the KAM theorem be ε 0 and δ be Under the above assumptions, with a probability of at least (1 − δ), a set of invariant tori exists for the trained model H NN .
Proof. It is confirmed by a straightforward calculation that if δ is given as described above, it holds that sup |H(u) − H NN (u)| < ε 0 , and hence the assumption of the KAM theorem is satisfied.
Remark 6. As mentioned in Remark 1, the KAM theorem also shows that the invariant tori become larger when the perturbation becomes smaller. Hence, if the generalization error is small enough, the size of the tori is expected to be large.
Note that general Hamiltonian systems, and hence general HNNs, are not quasi-periodic. Therefore, a model that approximates a quasi-periodic Hamilton equation may be (in some sense) approximately quasi-periodic, but it is not necessarily strictly quasi-periodic. This theorem states that the trained model can be strictly quasi-periodic even if the training loss does not completely vanish.
Numerical Example: Learning the Zabusky and Kruskal Experiment As a numerical experiment, we trained a HNN 1 so that the dynamics of the KdV equation is learned by using the data from the experiment by Zabusky and Kruskal (1965), in which a nontrivial recurrence of initial states is reported.
The KdV equation is derived from the following energy function H: In fact, under the periodic boundary condition, the variational derivative is Then, the time evolution is expressed as a Hamiltonian equation: For spatial discretization, we used the forward and backward difference operators, respectively. The central difference operator D is their mean, specifically D = 1 2 (D f +D b ) and that for the second derivative is D 2 = D f D b = D b D f . Using these difference operators, the equation is semi-discretized as Following Zabusky and Kruskal (1965), we set the parameters to α = −1.0 and β = −0.022 2 , set the width of phase space to 2.0, and used the initial condition u(0, x) to u(0, x) = cos(xπ). We discretized the system with the spatial and temporal mesh sizes of ∆x = 0.1 and ∆t = 0.01. We obtained an orbit for 200 time steps from the initial condition using the fifth-order Dormand-Prince method with the absolute and relative tolerances of 10 −10 and 10 −8 .
We performed the experiments on an NVIDIA TITAN V with double precision. We employed a three-layered convolutional neural network with kernel sizes of 3, 1, and 1. The number of hidden channels was 200, the number of output channels was 1, the activation function was the tanh function, and each weight parameter was initialized as a random orthogonal matrix. We summed up the output in the spatial direction and obtained the global energy. We used the whole orbit at every iteration, and minimized the mean squared error of the time derivative as the loss function using the Adam optimizer with a learning rate of 10 −3 for 10,000 iterations; the error reached a maximum of 1.37 × 10 −3 . Given the true dynamics u, the absolute error between the energy function H and the neural network H NN was 1.31 × 10 −4 on average and 2.51 × 10 −4 at most.
Using the true model and the trained neural network, we also obtained orbits for 1100 time steps from the same initial condition, as shown in the second and third panels of Fig. 4, respectively. In the top panels, blue and orange lines denote the true state u and the state predicted by the trained neural network u NN at t =0.0, 2.0, 4.8, and 9.8. The bottom panel shows the energy function H given the predicted states u and u NN . Due to the non-zero training error, more waves incur a larger error. Nonetheless, at around t = 9.8, the true model and learned neural network reproduce sin waves, which are given as the initial condition, and the energy error is restored to zero; they exhibit quasi-periodic behaviors.

Concluding Remarks
We analyzed the behavior of HNNs with non-zero learning errors by combining the KAM theory and statistical machine learning. We investigated the approximation properties of deep energy-based models, including HNNs. More precisely, we proved the persistence of the quasi-periodic behaviors of integrable Hamiltonian systems with a high probability even when the loss function is not perfectly zero. Further, we provided a generalization error bound and universal approximation theorems for HNNs to ensure that the loss function can be sufficiently small for application of the KAM theorem. Meanwhile, in the recent research on this type of model, numerically integrated gradients are often used for training. Similar results should be obtained for such cases; however, rigorous discussion is needed.

Energy-Based Physical Models
Equation (1) [and its coordinate transformation (3)] and the related model (2) [and its coordinate transformation (7)] describe not only Hamiltonian systems but also various types of other phenomena. In fact, as shown in the theorem below, models similar to (1), where G is a skew-symmetric or negative-semidefinite matrix, have not only the energy conservation property but also the energy-dissipation property, depending on the matrix G. The same result for the models similar to (3) is obtained in a straightforward way. Theorem 11. Equation (12) admits the energy conservation law, if G is skew-symmetric, and it admits the energy dissipation law, Proof of Theorem 11. We use the chain rule to obtain Hence, when the matrix G is skew-symmetric, the energy conservation law holds: dH dt = 0 because v Gv = 0 holds for all vectors v when G is skewsymmetric. When the matrix G is negative semidefinite, the energy function is monotonically non-increasing: As mentioned above, equations of this form are used not only in Hamiltonian mechanics but also in many fields of mathematical modeling, for example, in the Landau theory and the phase-field method. Phenomena such as phase separation, crystal growth, and crack propagation are modeled using these theories (e.g., Caginalp (1986); Miehe, Hofacker, and Welschinger (2010); Steinbach (2009) ;Wheeler, Murray, and Schaefer (1993)). The above model also includes the semi-discretized partial differential equations. For example, the KdV equation, which is a model of shallowwater waves, where α, β are parameters. This equation can be written as where H(u, ∂u/∂x) is the Hamiltonian and δH/δu is the variational derivative of H, which is defined by where u x denotes ∂u/∂x. Equation (15) is an example of a Hamiltonian partial differential equation, in which ∂/∂x in front of δH/δu plays the role of the skew-symmetric matrix G. In fact, if ∂/∂x is discretized by using a central difference operator, ∂ ∂x Because the difference matrix D is skew-symmetric, the semi-discretized equation has the form of equation (7) with a skew-symmetric matrix G. Similarly, an equation with the form is semi-discretized to equation (7) with a negativesemidefinite G. In fact, if the ∂ 2 /∂x 2 in front of δH/δu is discretized, which is known to be negative semidefinite (Furihata and Matsuo 2010). The universal approximatio-genn properties and the generalization error analysis are also applied to this model in a straightforward way under the assumption that S is nondegenerate. In fact, the properties of the matrix S used in the theoretical analyses are (1) S is non-degenerate and (2) the norm of S is bounded. If S is a constant matrix, (2) is automatically satisfied and hence (1) is the only assumption of the matrix S. In a part of the theoretical analysis, we consider the errors in the gradient ∇H of the energy function. However, typically, the loss function of S∇H is minimized in the training process, and if S is degenerate, errors in ∇H may not be small even when those in S∇H are small. Hence, S should be non-degenerate.

Proofs of Other Theorems
Proof of Theorem 5 Proof. By substituting the equation, we obtain because S is skew-symmetric and hence for any vector v, v Sv = 0.
Proof of Theorem 6 We prove the following theorem, which is a generalization of Theorem 6.
Theorem 14. Let H : R N → R be an energy function for the equation and S is an N × N non-degenerate matrix. Suppose that the phase space K of this system is compact and the right-hand side ∂H/∂u is Lipschitz continuous. Suppose also that u is a C 1 -diffeomorphism. If the functions σ = 0 and ρ = 0 belong to S 1 2 (R), then for any ε > 0 there exist neural networks H NN with the activation functions σ and u NN with ρ for which Proof. We need to prove the approximation property for (∂u/∂x) −1 . From the assumption that ρ = 0 is in S 1 2 (R), there exists a function u NN that approximates ∂u/∂x. Because the determinant function of matrices is continuous, it is deduced that det ∂u NN /∂x = 0 and hence (∂u NN /∂x) −1 exists. Because the matrix inverse is also continuous, (∂u NN /∂x) −1 is also approximated by u NN .

Details on Numerical Experiments
Energy-based physical model with the coordinate transformations In this study, we considered a model with coordinate transformations: For a better understanding of the significance of the theoretical results presented herein, we performed numerical tests using this model.
Because the generalized momenta are difficult to observe, we suppose that the values of the state variables θ 1 , θ 2 and their derivatives φ 1 , φ 2 are given as data. In such a situation, HNNs (Greydanus, Dzamba, and Yosinski 2019) are not applicable, because this model can be used only when p 1 and p 2 are generalized momenta; in the case considered here, they are supposed to be unknown.
To confirm this, we tested the model where v 1 , v 2 are not the generalized momenta but the velocities v 1 =q 1 , v 2 =q 2 . The experiments were implemented using Python 3.8.5 with the packages PyTorch 1.7.1, NumPy 1.19.5, SciPy 1.6.0, autograd 1.3, and torchdiffeq 0.2.1. For data, we used numerical solutions to the equation of motion with the parameters l 1 = l 2 = 1.0, m 1 = 1, m 2 = 2, and g = 9.8 solved by SciPy odeint with 100 initial conditions randomly generated from the standard normal distribution. We numerically integrated each orbit on the time interval [0, 5], in which the computed values are evaluated at 100 points with an identical sampling rate. Then, the target data d dt are obtained by substituting the computed state variables into the right-hand side of the equation. The experiments were performed on an NVIDIA GeForce 2080 Ti. The energy function was modeled by using a two-layered fully connected neural network with 50 hidden units and tanh as the activation function. We minimized the mean squared error of the time derivative as the loss function using the Adam optimizer with a learning rate of 0.001 and a batch-size of 200 for 10,000 iterations.
Remark 7. For the coordinate transformations in the model, invertible neural networks should be used; however, we used a simple neural network here because this model is assumed in the theorem in the main text. In the theorem, we needed neural networks that can approximate given functions and their derivatives. Although the universal approximation property of invertible neural networks was recently proved (Teshima et al. 2020), this theorem shows the universal approximation property of the function values only, not of the derivatives, which are needed for computation of the Jacobi matrix in our model.
Examples of the predicted orbits are shown in Figure 3. The training losses were 13.6 for the HNN and 0.280 for the model with coordinate transformations. As shown in this figure, the naive model failed to capture the dynamics correctly. This is because the dynamics of θ 1 , θ 2 , φ 1 , andφ 2 cannot be described by equation (17). This illustrates that, to model physical phenomena by using a model of the form the choice of the coordinate system is important.
Mass-spring system. Due to the well-known chaotic behaviors of the double pendulum, the results, in particular the values of the losses, of the previous experiments are to a certain extent unstable, except for the Hamilton neural networks always failing. Therefore, second, we investigated the models in more detail using the simple mass-spring system depicted in Figure 5. The two mass points m 1 and m 2 are connected by springs, which respectively have the spring Figure 5: Target mass-spring system.
constants k 1 and k 2 and natural lengths l 1 and l 2 . This system is a Hamiltonian system with the energy function where q 1 , q 2 are the positions of the mass points and p 1 , p 2 are the momenta, which are defined by p 1 = m 1 v 1 , p 2 = m 2 v 2 , v 1 = dp 1 /dt, andv 2 = dp 2 /dt. Suppose that we do not know the exact values of m 1 and m 2 and only the positions q 1 and q 2 and their derivatives can be observed.
Although m 1 and m 2 may be estimated from the data, for evaluation of the models, we tried to model the dynamics only using q 1 , q 2 , and their derivatives.
We examined the naive model, the model with the coordinate transformation, and the neural ordinary differential equation (ODE, Chen et al. (2018) The experimental conditions were almost the same as those in the previous experiment, except for the batch size, which we set to 100. For data, we used numerical solutions to d dt which is equivalent to equation (18).
Examples of the predicted orbits are shown in Figure 6. While the model with coordinate transformations produced an orbit that was nearly identical to the ground truth, the naive model failed to predict the states. In fact, it should be impossible to rewrite equation (19) as (17) with a certain energy function; for example, when the system has just one mass point and the equation of motion is given by d dt , this can be transformed into a Hamiltonian system d dt with k 1 = k 1 /m 1 , m 1 and the mass m 1 = 1. Hence, for this system, the HNNs are applicable without knowledge of m 1 . However, for system (19), such a transformation cannot be applied. Meanwhile, in the model with the coordinate transformation, a coordinate transformation that makes the equation Hamiltonian is explored. The result obtained by the neural ODE is better than that by the naive model; however, the prediction error became larger as time increased. This is due to the non-existence of the energy conservation law for the neural ODE. In fact, the value of the energy function of the neural ODE model was increasing, as shown in Figure 7.
The losses for the above models are listed in Table 1 for 12 experiments for each model. The values of the loss functions of the neural ODE were very small and those of the model with coordinate transformations were not stable. Actually, the histories of the loss functions for this model were not monotonic at all. This is due to the non-uniqueness of the model with coordinate transformations; if a model is fit to the given data quite well, then so are its canonical transformations. Due to this feature of the model, after the loss function decreased, it sometimes increased significantly in the search for different coordinate transformations. The failure in the learning process can be detected as the large train-ing loss, and if the learning process is successful, the test loss is competitive with that of the neural ODEs. Moreover, the neural ODEs exhibit energy drift, as seen above. Therefore, the model with coordinate transformations is the most suitable for a long-term physics simulation in an unknown coordinate system.