HHD-GP: Incorporating Helmholtz-Hodge Decomposition into Gaussian Processes for Learning Dynamical Systems

Authors: hao  xu
jia  pan
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 18:00:56

Abstract: Machine learning models provide alternatives for efficiently recognizing complex patterns from data, but the main concern in applying them to modeling physical systems stems from their physics-agnostic design, leading to learning methods that lack interpretability, robustness, and data efficiency. This paper mitigates this concern by incorporating the Helmholtz-Hodge decomposition into a Gaussian process model, leading to a versatile framework that simultaneously learns the curl-free and divergence-free components of a dynamical system. Learning a predictive model in this form facilitates the exploitation of symmetry priors. In addition to improving predictive power, these priors make the model indentifiable, thus the identified features can be linked to comprehensible scientific properties of the system. We show that compared to baseline models, our model achieves better predictive performance on several benchmark dynamical systems while allowing physically meaningful decomposition of the systems from noisy and sparse data.

Main paper: 
 introduction: 
 A dynamical system describes how the state of a system evolves over time [64] . Data-driven modeling of dynamical systems has become a fundamental task in many modern science and engineering applications, such as physical emulation [21] and robotics control [13] . Mathematically, a dynamical system is often a set of first-order ordinary differential equations (ODEs) or, equivalently, a smooth vector field on a manifold [28] . Given the functional form of ODEs, classical data-driven methods typically involve optimizing their parameters [52] . However, for many complex systems it is practically difficult to determine the form of the ODEs governing their dynamics.Recent advances in machine learning (ML) focus on the use of neural networks [11] and nonparametric Bayesian models [63, 27, 26] for the black-box approximation of vector fields. Despite the rich expressive power, the physics-agnostic design of these models hinders their performance when applied to physical systems. To address this issue, a popular approach is to develop ML models that incorporate strong physical priors as inductive biases. Such prior knowledge commonly stems from basic physical principles related to certain differential invariants of vector fields. For example, in scenarios of learning Hamiltonian systems [23, 67, 54] and incompressible fluid dynamics [73, 33] , ML models are constructed to learn divergence-free (div-free) vector fields, as a consequence of conservation laws of energy or mass. By constraining the solution space, the powerful physical principles effectively improve the extrapolation performance of ML models and enhance their interpretability. However, enforcing the model's behavior to adhere to centain fundamental physical principles also restricts the applicability of the model. For example, a div-free vector field fails to describe a dynamical system with dissipation, but real-world dynamical systems always suffer from non-negligible dissipation.To develop a predictive model covering more dynamical systems, we explore supplementing the divfree vector field with a curl-free vector field. This is inspired by the Helmholtz-Hodge decomposition (HHD) [3, 44, 8] , which states that any sufficiently smooth vector field can be expressed as the sum of a curl-free vector field and a div-free vector field. HHD is widely used in the study of Navier-Stokes equations [45, 20, 9] , but in this work we explore its connections with more general dynamical systems. For example, as shown in Fig. 1 , HHD can be used to characterize the dynamics of a dissipative Hamiltonian system, where the div-free component represents its conservative dynamics and the curl-free component describes the friction-induced dissipation. In addition to dissipative Hamiltonian systems, HHD is also widely used to study chaotic systems and ocean dynamics. The connection between HHD and these dynamical systems is described in more detail in Appendix A. The vector field in phase space of a massspring system and its HHD.In this work, we construct Gaussian process (GP) priors for div-free and curl-free dynamics separately, resulting in an additive GP model in the form of HHD for learning vector fields of dynamical systems.Our resulting HHD-GP allows one to leverage more prior knowledge than modeling the system as a whole. In particular, we investigate its potential in exploiting priors of symmetries, motivated by the observation that the div-free and curl-free components of a dynamical system usually exhibit more symmetry properties than the system itself. For example, the damped mass-spring system in Fig. 1 (a) exhibits odd symmetry, but its div-free (Fig. 1 (b)) and curl-free (Fig. 1(c )) components additionally present rotation and translation symmetry, respectively. Therefore, we further build symmetry-preserving div-free GPs and curl-free GPs by exploiting closure of GPs under linear transformation. The symmetry prior not only improves the predictive performance of HHD-GP, but also makes it identifiable, thus identified div-free and curl-free features can be physically-meaningful. In particular, the learned div-free features are closely related to the energy of dynamical systems.The main contributions of this work are summarized as:• We introduce a GP prior (called HHD-GP) for learning the Helmholtz-Hodge decomposition of a general dynamical system.• We construct a symmetry-preserving extension of the HHD-GP that can learn physically meaningful representations of dynamical systems.• Experiments on several benchmark systems show that our model can both accurately predict their dynamics and HHDs from sparse, noisy observations. 
related work: 
 Learning with div/curl-free constraints Div-free vector fields are a focal point of mathematical physics and have been well exploited by machine learning models to learn system dynamics, with the most well-known examples being neural networks (NNs) [23, 67] and GPs [54, 59] for learning Hamiltonian dynamics. All Hamiltonian vector fields are div-free, but not vice versa. To learn more general div-free vector fields as solutions of the continuity equation, [57] introduced an NN architecture to parameterize a universal representation of div-free vector fields. Based on the same representation, we constructed div-free kernels for GPs from matrix-valued kernels, which is an extension of the method of constructing div-free kernels from scalar kernels [50, 41, 75] . According to Maxwell's equations, div-free kernels were combined with curl-free kernels in [72, 71] to model magnetic fields. They assumed direct access to noisy observations of the div-free and curl-free components separately, whereas in this work, our goal is to recover the individual components from noisy observations of their sum. A similar prediction problem was studied in [6] , they combined 2D div-free and curl-free GP priors in the form of HHD to reconstruct planar ocean current fields, and recovered their divergence. Their model has the same formulation as ours when the dimension of HHD-GP is two, but we further developed a symmetry-preserving extension of HHD-GP to solve its non-identifiability problem. Another similar work is Dissipative Hamiltonian neural network (D-HNN) [22] , which compensated HNN [23] with a curl-free part to model both conservative and dissipative dynamics simultaneously, but D-HNN is only applicable to dissipative Hamiltonian systems due to its construction.Learning with symmetry Symmetries are another important aspect of priors that can be incorporated into machine learning models. Motivated by the success of the translation-invariant NNs [38] , network architectures with symmetries to more general transformations have been proposed, such as steerable CNNs [74, 10] and graph NNs [46, 61] equivariant to Euclidean symmetries. They achieved great success in improving generalization of the models. With the same motivation, kernel methods incorporating symmetries have also been developed. To make predictions invariant to transformations of inputs, [24] constructed kernels using Haar integration. And based on similar integration technique, [56] developed the Group Integration Matrix Kernel (GIM-kernel) to learn equivariant functions, which was later used by [58] to learn dynamical systems with symmetries. In this work, we constructed GP models to impose Euclidean symmetries to div/curl-free vector fields, which to the best of our knowledge has not been explored by the machine learning community.3 BackgroundThe group of all rotations in R n is represented by the set of special orthogonal matricesSO (n) = {R ∈ O (n) | det R = 1}, where a rotation matrix R transforms a point by L R (x) = Rx. 3. Reflection: Reflections in R n forms subgroups of the orthogonal group O (n). Reflections correspond to mirror symmetries. They mirror points across a hyperplane. For a hyperplane with a unit normal vector n, the action of a reflection is defined asL n (x) = x -2(x • n)n. 
hhd and problem setup: 
 We consider an autonomous system governed by the following ODEs:EQUATIONwhich defines a vector field by assigning a vector f (x) ∈ R n to every state x ∈ R n . We assume that the vector field f ∈ L 2 (R n , R n ) is smooth, and any such vector field can be decomposed into the sum of a curl-free vector field f curl : R n → R n (∇ ∧ f curl = 0, ∀x ∈ R n ) and a divergence-free vector field f div : R n → R n (∇ • f div = 0, ∀x ∈ R n ), according to the Helmholtz-Hodge decomposition (HHD) [44, 8] . In this work, we are interested in learning f , f curl and f div simultaneously from a collection of noisy observations denoted by D = {(x i , y i )} m i=1 , with a noisy observation y i at a state x i given byEQUATIONwhere an additive noise ϵ ∈ R n follows a zero-mean Gaussian distribution defined by a covariance matrix Ω = diag σ 2 1 , . . . , σ 2 n ∈ R n×n modeling noise variance in each output dimension.We are interested in using Gaussian processes (GPs) to infer unknown vector fields. A GP is a stochastic process commonly used as a distribution for functions, assuming that any finite number of function values has a joint Gaussian distribution [53] . To learn an unknown vector field f : R n → R n , we assume a vector-valued GP prior as f (x) ∼ GP (0, κ (x, x ′ )), where the mean of the function values is set to zero, and the covariance is captured by a matrix-valued kernel κ : R n × R n → R n×n , whose (i, j)-th entry expresses the correlation between the i-th dimension of f (x) and the j-th dimension of f (x ′ ). In the GP framework, the kernel controls the properties of possible functions under a GP, leading to various efforts of problem-specific design of kernels [76, 16, 18] . And to be a valid covariance function, the kernel should be symmetric and positive semidefinite [78] .GPs provide a Bayesian non-parametric approach for solving regression tasks. According to the GP prior, function values at inputs X = [x 1 , . . . , x m ] T are jointly distributed as N (f (X) ; 0, K), T given by Eq. 2 can be calculated bywhere K = [κ (x i , x j ) ∈ R n×n ] m i,j=1EQUATIONwhere K * = K + Σ, Σ = Ω ⊗ I m is a diagonal matrix whose elements are the variance of the observation noise. Training a GP model refers to maximizing the log of Eq. 3 to optimize the kernel parameters and the noise variance. Then by conditioning on these observations using Bayes' rule, the predictive posterior for a new state x * is still Gaussian with its mean µ and variance v given byEQUATIONwhere k = [κ (x 1 , x * ) , . . . , κ (x m , x * )] T ∈ R mn×n . The derived mean function is used for regression results, and the associated variance quantifies prediction uncertainty. Due to their nonparametric nature, GPs can self-adapt to the complexity of the target function based on the data provided, without being restricted to specific parametric forms.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
hhd-gp model: 
 We consider the problem of learning a continuous-time dynamical model in the form of HHD (Eq. 1) with GPs. HHD points out the prevalence of an additive structure in dynamical systems, so the key idea here is to exploit two GPs to model f curl and f div , respectively, f curl ∼ GP (0, κ curl ), f div ∼ GP (0, κ div ). Then the sum of these two GPs results in a new GP modeling the dynamical system f , with a new kernel function defined as the sum of the curl-free and divergence-free ones. And the additive kernel κ = κ curl + κ div inherits the symmetric and positive-semidefinite properties of κ curl and κ div , so the GP predictor (Eq. 4) is valid for the additive GP model. And the additivity of the kernels implies the additivity of GP means, so the mean function µ in Eq. 4 can be split into a curl-free part µ curl and a divergence-free part µ div , and we haveEQUATIONwhere K * = K curl + K div + Σ. It can be seen that the effects of f curl and f div can be treated as observation noises for each other, so their prediction variances at x * are obtained byv curl (x * ) = κ curl (x * , x * ) -k T curl K -1 * k curl , v div (x * ) = κ div (x * , x * ) -k T div K -1 * k div . (6)Consequently, observations of f (x) given by Eq. 2 can be used to make predictions for its hidden components f curl and f div . So we then construct GPs with realizations in the space of curl-free and div-free vector fields.In the following parts, we construct the kernels κ curl and κ div from the representations of f curl and f div respectively, exploiting closure of GPs under linear transformation. Let R = R [∂ x1 , . . . , ∂ xn ] be the polynomial ring in the partial derivatives, and L ∈ R a×b be a matrix of differential operators acting on functions g : R n → R b distributed as GP (µ (x) , κ (x, x ′ )). Then, the transformation of g under L is again distributed as a GP withEQUATIONwhere L x and L x ′ denote the operation of L on the first and second argument of κ (x, x ′ ), respectively [29, 1, 36, 7] . To make L x κ (x, x ′ ) L T x ′ a valid covariance function, its underlying kernel κ (x, x ′ ) must be twice differentiable in R n . See Appendix B for details on linear operations on GPs.The gradient operator(∇ := [∂ x1 , . . . , ∂ xn ] T ∈ R n×1) defines a surjective mapping from the space of smooth scalar fields to the space of curl-free vector fields [14] , so f curl can be represented byf curl = ∇V , where V ∈ C ∞ (R n , R) is called the scalar potential of f curl . Since the gradient operation defines a linear transformation, if a GP with a scalar kernel κ V is assumed on V , the distribution of f curl is again a GP. According to Eq. 7, the curl-free GP over f curl is given byEQUATIONwhere κ curl : R n × R n → R n×n is a matrix-valued kernel constructed by the Hessian of the scalar kernel κ V , consisting of all second-order partial derivatives in x and x ′ , with the entry in the i-th row and j-th column given by[κ curl (x, x ′ )] i,j = Cov ∂V (x) /∂x i , ∂V (x ′ ) /∂x ′ j = ∂ 2 κ V (x, x ′ ) /∂x i ∂x ′ j .By this construction, if κ V induces a GP with realizations dense in C ∞ (R n , R), the set of realizations of GP (0, κ curl ) is dense in the space of curl-free vector fields, because a surjective mapping maps dense sets to dense sets.A div-free vector field can be constructed from a skew-symmetric matrix field [5, 32, 57] . Specifically, let A : R n → R n×n be a skew-symmetric matrix-valued function, then a div-free vector field f div can be represented by taking row-wise divergence of A, i.e.,EQUATIONwhere A i : R n → R n is the i-th row of A. The skew-symmetric matrix field A of size n × n can be compactly represented by m = n (n -1) /22 scalar functionsEQUATIONwhere Φ ij ∈ R n×n is a matrix with its (i, j)-th entry equal to 1, (j, i)-th entry equal to -1, and all other entries equal to 0. Then, the div-free vector field given by Eq. 9 can be reformulated as a linear transformation:EQUATIONwhereψ ij = Φ ij ∇ ∈ R n×1is a column vector obtained by linearly transforming the gradient operator. The m column vectors ψ ij are aggregated in a matrix Ψ ∈ R n×m , and the m corresponding scalar functions u ij are collected in a vector-valued function u : R n → R m . Ψ [•] is a matrix of linear differential operators, so to use a GP to model f div , we can proceed by assuming a GP prior over u ∼ GP (0, κ u : R n × R n → R m×m ), then based on the closure of GPs under linear transformation (Eq. 7), the GP prior over f div (Eq. 11) is constructed byEQUATIONwhere κ u is a scalar-valued kernel for two dimensional systems (n = 2, m = 1), and is a matrixvalued kernel for n > 2. Notice that this div-free kernel can be treated as a generalization of the div-free kernel derived from a scalar kernel [50, 41, 75] , which can be recovered by setting κ u = κ • I. Therefore, it is more expressive and flexible. Theoretically, the GP model given by Eq. 12 can be used to approximate arbitrary div-free vector fields, because the representation given by Eq. 9 is shown to be maximally expressive (i.e. universal) in [57] .With the curl-free and div-free kernels, our objective is to learn physically interpretable representations of a dynamical system based on the HHD-GP model. However, the HHD is always not unique due to the existence of harmonic components f harm (vector fields satisfying both ∇ ∧ f harm = 0 and ∇ • f harm = 0, e.g., constant vector fields). For the HHD of a dynamical system with the true functional decompositionEQUATIONis a valid HHD for arbitrary f harm , which thus makes the HHD-GP model non-identifiable, meaning that from the same training data, we may learn different decompositions giving the same predictions. This is not desirable because we expect the learned dynamical model to be interpretable: the curl-free and div-free components f curl , f div are physically meaningful.To mitigate the identifiability problem in additive regression models, an effective method is to impose constraints on their component models [16, 17, 47, 42] . The imposed constraints can affect the decomposition results of the additive models. Therefore, to ensure that HHD-GP can produce a scientific decomposition, we desire constraints that respect the inherent characteristics of dynamical systems. And, as another primary goal, incorporating prior knowledge of a system into a GP model can also improve its prediction accuracy and learning efficiency. Therefore, in the next section, we present how to impose symmetry-based constraints on the curl-free and div-free GP models.5 Symmetry Constraints 
equivariance and invariance: 
 Symmetry is a fundamental geometric property prevalent in dynamical systems in natural [39] , and is usually described by the concept of equivariance and invariance: Definition 5.1 (Equivariance and Invariance). Let G be a group acting on R n through a smooth mapEQUATIONwhere L g (x) := L (g, x), and J Lg denotes the Jacobian matrix of L g . Then, G is termed the symmetry group of the dynamical system. In particular, if J Lg is the identity matrix (i.e., f • L g = f , ∀g ∈ G), the dynamical system f is said to be G-invariant.From the equivariance condition (Eq. 14) of the vector field, it follows the system's trajectory commutes with the action map. For vector fields on R n , the symmetry group G is commonly a subgroup of the Euclidean group E (n), which comprises all intuitive geometric transformations in R n (see Appendix C for a brief introduction). The symmetry constraints refer to that we expect the learned curl-free and div-free vector fields to be G-equivariant. With the representation of their GP models, we demonstrate that the symmetries can be enforced by designing suitable kernels.The curl-free GP (Eq. 8) is constructed by transforming another GP over a potential function, implying that we can impose constraints of symmetry on the curl-free GP by designing a suitable potential GP. Therefore, we start by exploring how to construct potential functions to obtain curl-free vector fields with the desired equivariance. As expected, the following theorem holds:Theorem 5.2. Let G be a Euclidean group or its subgroup, and let V : R n → R be a G-invariant scalar function. Then, the curl-free vector field f curl : R n → R n defined byf curl (x) = ∇V (x) is G-equivariant.See Appendix D.1 for the proof. Theorem 5.2 shows that a G-invariant scalar potential V can yield a G-equivariant gradient field, indicating that if any realization V of GP (0, κ V ) is constrained to be G-invariant, then its pushforward GP over ∇V ∼ GP 0, ∇ x κ V ∇ T x ′ can induce the space of G-equivariant curl-free vector fields.It is obvious that a G-invariant scalar potential V can be constructed by integrating some non-invariant function h : R n → R over the symmetry group: V = G (h • L g ) dg, where the measure dg is called Haar measure, which exists for locally compact topological groups and finite groups. Therefore, by assuming that h is distributed as h ∼ GP (0, κ h ), we can construct the GP prior over the G-invariant scalar potential as V ∼ GP (0, κ V ), with its kernel κ V given byEQUATION)This kernel is called the Haar-integration kernel [24] . While it provides a general method for constructing kernels for G-invariant functions, the double integral can be computationally expensive.If the kernel κ h is invariant to any g ∈ G in the sense that κ (x, x ′ ) = κ (L g (x) , L g (x ′ )) 3 , a complexity reduction of Eq. 15 by one square-root can be performed byEQUATIONwhere |G| = G dg, and it denotes the cardinality of G when the group is finite.To incorporate the equivariance condition (Eq. 14) into realizations of the div-free GP (Eq. 12), we construct the skew-symmetric matrix field A from a vector-valued function. Specifically, given a smooth vector fieldv ∈ C ∞ (R n , R n ), A is constructed by A = J v -J T v, where J v denotes the Jacobian of v with its (i, j)-th entry given by ∂v i /∂x j . Then the component function u ij in Eq. 10 is given by u ij = ∂v i /∂x j -∂v j /∂x i . By this construction, the symmetry of the div-free vector field f div is governed by the symmetry of its vector potential v. In particular, a G-equivariant v can produce a G-equivariant f div , and it is formalized in the following theorem: Theorem 5.3. Let G be a Euclidean group or its subgroup, and let v : R n → R n be a G-equivariant vector field. Then the divergence-free vector fieldf div = [∇ • A 1 , . . . , ∇ • A n ] T is G-equivariant,where A i denotes the i-th row of the skew-symmetric matrix-valued functionA = J v -J T v .See Appendix D.2 for the proof. By this theorem, we then proceed by assuming a GP prior over the vector potential v ∼ GP (0, κ v ), and to constrain v to be G-equivariant, we build its kernel κ v ∈ R n×n in the form of the Group Integration Matrix kernel (GIM-kernel) [56, 55] , which is constructed by:EQUATION)where κ is some arbitrary scalar-valued kernel satisfying κ (x,x ′ ) = κ (L g (x) , L g (x ′)) for all g ∈ G. The GIM-kernel spans a Reproducing Kernel Hilbert Space (RKHS) of functions with the desired equivariance [55] . So we can then use κ v (Eq. 17) to construct the GP prior over u in Eq. 11, where the covariance between components u ij and u kq is given byEQUATIONFinally, this matrix-valued kernel κ u is transformed by Eq. 12 to construct the div-free GP, of which the realizations are guaranteed to be G-equivariant div-free vector fields, according to Theorem 5.3. 
experiments: 
 We evaluated the performance of our proposed method in several representative physical systems.Through the experiments, we found that our model can not only accurately capture the system dynamics, but also predict correct decompositions.We first evaluated our method on a damped mass-spring system and a damped pendulum. Their governing equations are detailed in Appendix E.1. We generated the training data {(x, ẋ)} by randomly sampling states x in their phase space, and each of their derivative observations ẋ is corrupted by an additive Gaussian noise with a standard deviation of 0.05.The models are first evaluated in terms of learning ODEs. Specifically, our evaluation focused on the accuracy of the models in predicting state derivatives, as measured by the root mean squared error (RMSE), and their ability to accurately predict state trajectories over time, as indicated by the valid prediction time (VPT). In addition to these metrics for evaluating the regression results, we further use the mean negative log likelihood (MNLL) to evaluate the prediction uncertainty provided by the GP models. These evaluation metrics and the generation of test data are detailed in Appendix F. We compared our models, HHD-GP and its symmetry-preserving extension, SPHHD-GP, with Dissipative Hamiltonian neural network (D-HNN) [22] , GPs involving div-free kernels for learning Hamiltonian dynamics (Div-GP) [54, 59] , and GPs with Group Integration Matrix Kernels (GIM-GP) [56] that can incorporate symmetries. Another baseline is GPs with independent kernels (Ind-GP), which model each dimension of a dynamical system with an independent scalar GP. See Appendix G for the implementation details of these models. The results for 20 training data are shown in Table 1 , and the results for an increasing number of training data are provided in Appendix J.1. The performance of Div-GP is limited because it can only model conservative dynamics. HHD-GP improves its performance by compensating with a curl-free kernel, which offsets the strong inductive bias imposed by the div-free kernel. And the performance of HHD-GP is better than that of another HHD-based model, D-HNN, because the low data efficiency of NNs makes it hard for D-HNN to capture dynamics using noisy and sparse training data, so actually the performance of D-HNN is worse than either of the GP methods. As another model without inductive bias, Ind-GP performs similarly to HHD-GP in most cases. Then, by incorporating symmetry priors into GPs, GIM-GP performs better than the above models but not as well as SPHHD-GP, because learning in the form of HHD allows the model to exploit more implicit symmetries in the dynamical systems. SPHHD-GP performs best overall in learning ODEs. Appendix J.2 presents the plots of trajectory predictions for each system. Another advantage of our model is that it can decompose the dynamics into its div-free and curl-free components. According to Eq. 19, the div-free component can be used to recover the system's Hamiltonian, as long as the correct form of HHD is learned by our model. To show this, we evaluated our model by another task: we predicted the Hamiltonian Ĥ (x t ) (energy) along the system's trajectory {x t }. For HHD-GP and SPHHD-GP, the Hamiltonian is predicted from a joint GP prior over the Hamiltonian and its underlying dynamics (see Appendix G.7.1 for details). Fig. 2 shows the predicted energy evolution of the two systems, where we can find that D-HNN and HHD-GP fail to provide physically plausible results, because the energy should continue to decrease due to friction, but they provide oscillating predictions, along with significant variances. In contrast, predictions of SPHHD-GP are highly accurate and closely aligned with the true values. One reason is that the symmetry priors used by SPHHD-GP improves the generalization performance of the model, but more importantly, the priors solve the problem of non-identifiability suffered by HHD-GP and D-HNN. See Appendix J.3 for an visualization of predicted decompositions, which shows that although the models capture similar system dynamics, they can learn completely different decompositions. Consistency is a necessary condition for a learning model to be identifiable and refers to the property that its parameter estimates should converge to the true values as the amount of data increases [19, 40] . Therefore, to further explore the non-identifiability problem, we provide the RMSE of energy prediction with an increasing number of training data in Fig. 3 .The result shows that HHD-GP and D-HNN always generate significant errors in predicting energy evolution, meaning that the model parameters cannot converge as the amount of data increases. In contrast, the energy errors of SPHHD-GP exhibits a decreasing trend, which reflects the consistency of SPHHD-GP. In addition to this empirical validation, in Appendix H we provide a theoretical verification that for the dynamical systems in our experiments, the non-identifiability of our model (i.e. non-uniqueness of HHD) is solved through forced symmetries.Accurate prediction of energy demonstrates the interpretability of learned div-free features. To further show the interpretability of curl-free features, in Appendix J. 4 we added experiments of adapting the learned models to predict dynamics with unseen friction coefficients. And Appendix J.5 presents some additional experimental results of increasing the noise in training data, showing that SPHHD-GP maintains high gain in large noise compared to the baselines. Learning the correct HHD also has important implications for studying chaotic systems. [60] stated that the div-free component of a chaotic system is always orthogonal to the gradient of its energy function, i.e., ∇H T f div (x) = 0. This energy function can be used to synchronize two chaotic systems [60] , analyze their stability [77] , and design energymodulated controllers [43] . This experiment aims to learn the Chua circuit [49] , which is a chaotic system with applications in various fields. Its governing equation and the HHD are detailed in Appendix E.2. The models are trained on 100 randomly sampled data, corrupted with additive Gaussian noise (standard deviation: 0.05). The generation of test data is the same as in the experiment in Section 6.1. The results are given in Table 2 , which shows that SPHHD-GP consistently outperforms the other approaches. Appendix J.6 plots the trajectory and energy predictions. These results again confirm the advantage of incorporating symmetry constraints into our method. Please note that D-HNN is not applicable to this system because the div-free part of D-HNN only applies to Hamiltonian systems.To investigate the performance of our model in learning real-world dynamical systems, we then evaluated our model on an ocean current field, which is a complex dynamical system intricately governed by the interplay of multiple factors such as Earth's rotation, wind patterns, temperature gradients, and coastal interactions. Based on sparse observations of buoy velocities, oceanographers are interested in estimating ocean currents away from buoys and identifying divergences of the ocean current field.In this experiment, we used a dataset from [6] , containing 1183 velocity data points from 12 buoys distributed across the northern Gulf of Mexico, as shown by the red arrows in Fig. 4 . We compared our model with the Helmholtz-GP model [6] , which combines 2D div-free and curl-free kernels in the form of HHD to reconstruct planar ocean current fields, but without considering the severe non-identifiability problem caused by the non-uniqueness of HHD. To make the HHD of an ocean current field unique, an effective way in the field of fluid dynamics is to enforce the parallel boundary condition, which states that the HHD of a vector field in a bounded domain is unique if its div-free component is parallel to the domain boundary (cf. page 36 in [12] , Section 5.1 in [8] ). For our proposed model, the method of incorporating symmetry can be used to enforce the parallel boundary condition. Specifically, if a vector field has mirror symmetry with respect to a hyperplane, then the vectors on this hyperplane must be parallel to the hyperplane. Therefore, to make our model identifiable, the div-free kernel in our model was constructed to incorporate mirror symmetry with respect to the rectangular boundary of the ocean current field. To reduce the computational complexity, the kernels were fitted into a sparse GP framework [62] , which reduced the computational complexity from O m 3 to O M 2 m , where m = 1183 is the number of data points and M = 200 is the number of inducing points. See Appendix I for the details of computational complexity of our model. The predictions of the ocean current field and its HHD are shown in Fig. 4 . Although the ground truth is not available for this field of real-world ocean currents, we can still assess performance against oceanographers' expert knowledge (cf. Section I.3 in [6] ). For the predictions of ocean currents (the first column of Fig. 4 ), the oceanographers expect to see continuous currents with no sharp deviations. However, we can clearly observe that the prediction of Helmholtz-GP shows an abrupt drop in positions away from observed data, especially in the lower left, upper right and center regions. In contrast, our model (SPHHD-GP) presents a more continuous current prediction. The predicted divergence field (i.e., the scalar potential of the curl-free component) is provided in the third column of Fig. 4 . The oceanographers expect to find a rich structure in divergence predictions. However, the Helmholtz-GP fails to recover the divergence field and instead predicts an almost constant curl-free vector field, as shown in the third block of Fig. 4(a) . This is caused by the non-identifiability of the Helmholtz-GP, since harmonic components usually exist in the form of constant vector fields (cf. Section 4.3). Our model, instead, recovers the rich structure of the divergence field. Therefore, we can conclude that our model provides more realistic ocean current predictions and divergence identification, offering a better alternative for the simulation of ocean dynamics. 
conclusion and future work: 
 Our work develops an additive GP model whose component is either free of divergence or of curl, the two most ubiquitous differential invariants of vector fields in natural, and we constrain the div/curl-free kernels to preserve desired symmetries. These symmetry-preserving kernels not only improve the accuracy of predictions but also make the model identifiable, thus a physically meaningful decomposition can be predicted. Our future work is to extend our model to exploit the connection of HHD with more dynamical systems. For example, there are recent advances in using the HHD to construct Lyapunov functions [65] . So, our model has potential to achieve good performance in learning stable dynamics. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We used a single Nvidia GeForce GTX 3050 Ti GPU in the experiments, as mentioned in the first paragraph of Appendix G.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: [TODO]Guidelines: We have discussed the societal impacts in Appendix K.• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: The code packages and datasets are properly cited in the paper.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: This section provides a brief introduction to the connections between the Helmholtz-Hodge decomposition (HHD) and some dynamical systems needed to understand this work.We begin with a brief review of Hamiltonian mechanics [4] . For a dynamical system with N degrees of freedom, the Hamiltonian formalism describes the system by defining a scalar function H (x) known as the Hamiltonian, where the system state x = (q, p) ∈ R 2N is described by generalized coordinates q ∈ R N and p ∈ R N in the phase space, corresponding to generalized position and momentum, respectively. The time evolution of (q, p) is governed by the symplectic gradient of its Hamiltonian, i.e., q = ∂H ∂p , ṗ = -∂H ∂q . If the Hamiltonian is not explicitly time-dependent, [ q, ṗ]T defines a divergence-free vector field ( ∂ ∂q ∂H ∂p -∂ ∂p ∂H ∂q = 0) whose flows conserve the Hamiltonian. By interpreting the Hamiltonian as energy, the Hamiltonian vector field can model systems with energy conservation. However, real-life systems often suffer from energy dissipation. In physical systems governed by an autonomous ODE, any energy variation occurs along a volume change in phase space and vice versa. Therefore, to account for energy dissipation, a term responsible for the volume contraction of the phase space can be introduced into the Hamiltonian dynamics. This term represents the energy lost by the system due to various dissipative forces and is typically modeled by gradient fields of some scalar functions D (•), so the motion equations of the dissipative Hamiltonian system can be given byEQUATIONIn a dissipative Hamiltonian system, the Hamiltonian still represents the system's total energy, but the additional damping term causes the system to lose energy over time. One common example of the damping term is induced by the Rayleigh function, D (p) = -1 2 p T Qp, where Q ∈ R N ×N is a symmetric positive-definite matrix called the Rayleigh dissipation matrix. The Rayleigh function provides an elegant way to include dissipative forces-such as friction, air resistance, and viscosity-in the context of Hamiltonian mechanics. When the Hamiltonian H (•) represents mass instead of energy, Eq. 19 can also be used to model 2-dimensional compressible fluid dynamics, such as the dynamics of ocean current fields.The dissipative Hamiltonian system (Eq. 19) is in an explicit form of HHD, because the Hamiltonian vector field is divergence-free and the dissipative field ∇D is curl-free. However, as its divergencefree part is governed by the Hamiltonian equations, Eq. 19 can only be used to describe a subset of even-dimensional systems. To extend its scope to study chaotics systems, [60] proposed to approach energy for dimensionless dynamical systems by using HHD as follows,EQUATIONwhere f div and f curl are divergence-free and curl-free components of a dynamical system, respectively. According to this criterion (Eq. 20), HHDs of chaotic systems were approached to calculate the energy for synchronizing two chaotic systems [60] , to analyze the stability of chaotic systems [77] , and to design energy modulation-based controllers [43] .Gaussian processes (GPs) are closed under linear transformations (cf. Lemma 2.2 in [36] , Lemma 2.1 in [25] ). Let L be a linear operator acting on realizations of g ∼ GP (µ g (x) , κ g (x, x ′ )), then under the assumption that L commutes with expectation, the mean of Lg is given byEQUATIONand the covariance is such thatEQUATIONA common application of this technique is to construct GPs with realizations in the solution set of linear differential equations, assuming that L is a linear differential operator [29, 35, 36, 7, 25] .And similarly, we make use of the closure of GPs under linear differential operators to construct the curl-free and div-free kernels.Although the transformed kernel (Eq. 22) has been widely used, its validity as a covariance function is rarely discussed in existing works. In a GP framework, a kernel κ (x, x ′ ) : R n × R n → R m×m is a valid covariance function if it is:(i) Symmetric, i.e., κ (x, x ′ ) = κ (x ′ , x) T , ∀x, x ′ ∈ R n ,(ii) Positive semidefinite, i.e., ij c T i κ (x i , x j ) c j ≥ 0, for any finite set {x i } ⊂ R n and {c i } ⊂ R m . These two conditions for a kernel can be verified if and only if there exists a feature map ϕ (x) such that κ (x,EQUATIONTherefore, the transformed kernel L x κ g (x, x ′ ) L T x ′ is guaranteed to be a valid covariance function provided that its underlying kernel κ g (x, x ′ ) is. However, when L is a differential operator, κ g should be twice differentiable, which is satisfied by most of the standard kernels, such as the squared exponential kernel.This section gives the basic definitions about the Euclidean group. In the context of this work, the most important example of a symmetry group is the Euclidean group E (n) or its subgroups. The set of all elements in E (n) can be denoted asEQUATIONwhereO (n) = A ∈ R n×n | AA T = I is the orthogonal group. Any element g = (A, b) ∈ E (n)represents a translation followed by an orthogonal transformation, the action of g on a point x ∈ R n is given by a linear mapping:EQUATIONTherefore, E (n) comprises all isometries of a Euclidean space, i.e. for all x, x ′ ∈ R n and g ∈ E (n), we haveEQUATIONwhere ∥•∥ is the Euclidean norm. All intuitive geometric transformations in R n can be described by subgroups of E (n), such as 1. Translation: The group of all translations in R n is denoted by (R n , +). For any v ∈ R n , A translation is a transformation that moves a point x ∈ R n by L v (x) = x + v.This section restates and proves the theorems in Section 5, which give the theoretical foundations that we use to enforce symmetry constraints to our HHD-GP model.D.1 Proof of Theorem 5.2 Theorem 5.2. Let G be a Euclidean group or its subgroup, and let V : R n → R be a G-invariant scalar function. Then, the curl-free vector field f curl : R n → R n defined byf curl (x) = ∇V (x) is G-equivariant.Proof. To prove that ∇V (x) is G-equivariant, by Definition 5.1 we need to show that for any x ∈ R n and g ∈ G, we haveEQUATIONwhere L g (x) : R n → R n is the action of g on x, and J Lg is the Jacobian matrix of L g . To show this, we first recall the chain rule of the gradient operator,EQUATIONwhere J T Lg denotes the transpose ofEQUATIONNow taking the gradient of both sides, and combining with the chain rule, we obtainEQUATIONThen by multiplying J Lg on both sides, we haveEQUATIONConsidering G is a subgroup of the Euclidean group, the action of any group element g ∈ G on R n can be represented byEQUATIONwhere A ∈ O (n) is an orthogonal matrix. Therefore, J Lg J T Lg = AA T = I, where I is an identity matrix, so we have J Lg ∇V = ∇V • L g , which means that the curl-free vector field f curl = ∇V is G-equivariant, as desired. D.2 Proof of Theorem 5.3 Theorem 5.3. Let G be a Euclidean group or its subgroup, and let v : R n → R n be a G-equivariant vector field. Then the divergence-free vector fieldf div = [∇ • A 1 , . . . , ∇ • A n ] T is G-equivariant,where A i denotes the i-th row of the skew-symmetric matrix-valued function A = J v -J T v . Proof. To prove this theorem, we need to show that the divergence-free vector field f div satisfies the G-equivariance condition, i.e., for all g ∈ G and x ∈ R n ,EQUATIONwhere L g (x) : R n → R n is the action of g on x, and J Lg is the Jacobian matrix of L g .Define the divergence-free vector fieldf div = [∇ • A 1 , ..., ∇ • A n ] T, where A i is the ith row of the skew-symmetric matrix field A = J v -J T v , and J v is the Jacobian of some vector field v : R n → R n . Given that v is G-equivariant, for all g ∈ G and x ∈ R n we haveEQUATIONNow computing the Jacobian of both sides, by the chain rule of the Jacobian, we obtainEQUATIONThen applying J T Lg on both sides, we haveJ T Lg (J v • L g ) J Lg = J T Lg J Lg J v .Since G is a subgroup of the Euclidean group, the Jacobian of the action of any group element g ∈ G on R n is an orthogonal matrix, i.e., J Lg ∈ O (n). Therefore, J T Lg J Lg = I, where I is an identity matrix, so it holds thatEQUATIONwhere the equation at the (i, j)-th entries of both sides is given byEQUATIONThen by substituting A ij into the construction of f div , for any g ∈ G and x ∈ R n , we haveEQUATIONThen, applying J Lg on both sides, we obtain J Lg f div = f div • L g , which completes the proof.This section reviews the physical systems used in our experiments, listing their governing equations, HHDs and symmetry properties.The connection between HHD and dissipative Hamiltonian systems is reviewed in Appendix A.1. In our experiments, the models are evaluated with the following two dissipative Hamiltonian systems.Damped Mass Spring A damped mass-spring system is a mass attached to a spring that oscillates periodically around an equilibrium position. Its Hamiltonian in natural units is given byEQUATIONwhere q ∈ R is its position, and p is the momentum conjugate to q. The Hamiltonian represents the total energy of the oscillator. Without energy dissipation, the motion of the oscillator is described by the Hamiltonian vector field (div-free vector field f div ) derived from the Hamiltonian:EQUATION)Then we show that the Hamiltonian vector field has SO(2)-equivariance. Apply a rotation by an angle θ that transforms any state vector [q, p] T to [q ′ , p ′ ] T , where q ′ = q cos θ-p sin θ, p ′ = q sin θ+p cos θ, then for all θ ∈ R, the equivariance condition (Eq. 14) can be easily verified byEQUATIONTherefore, the Hamiltonian vector field is equivariant with respect to the 2D rotations. However, this rotation symmetry is broken when the system suffers from energy dissipation, i.e., a dissipative field is added to the motion equations. The dissipative vector field (curl-free vector field f curl ) is induced by the Rayleigh function D (p) = ρp 2 /2:EQUATION)where ρ is the friction coefficient and we set it at 0.1. This dissipative field has the symmetries of the q-axis translation (f curl (q + g, p) = f curl (q, p), ∀g ∈ R) and the rotation by an angle of π radians (the so-called odd symmetry, i.e., f curl (-q, -p) = -f curl (q, p)). Then by summing the curl-free vector field (Eq. 42) and the Hamiltonian vector field (Eq. 40), the dynamics of the damped harmonic oscillator is characterized byEQUATIONwhich only demonstrates the odd symmetry.Damped Pendulum A damped pendulum is a physical system consisting of a weight suspended from a pivot, subjected to a resistive force that gradually reduces its oscillation over time. By defining its state in the phase space through [q, p] T ∈ R 2 , the dynamics of a damped pendulum and its HHD can be given byEQUATIONwhere the Hamiltonian H and the Rayleigh function D are given byEQUATION)In our experiments, the gravitational constant g, the mass m, and the pendulum length l were set as g = 3 and m = l = 1, and the friction coefficient ρ was set as ρ = 0.1. The dynamics of the damped pendulum, as well as its div-free and curl-free components, exhibit the odd symmetry, but the curl-free part additionally exhibits the translation invariance along the q-axis.The Chua circuit [49] is a simple electronic circuit that exhibits chaotic behavior, and it has applications in various fields, such as secure communication systems and random number generators. The ODE of a Chua circuit with its HHD is given byEQUATIONwhere x, y, z are the phase space variables, and α, c, β are the system parameters, which were set to 10, -0.143, and 16 respectively in the experiments. The Chua circuit is equivariant under a π rotation about the origin: L π : (x, y, z) → (-x, -y, -z). If the Chua circuit is decomposed by the HHD in Eq. 46, more symmetries are exhibited. Specifically, its div-free component inherits the odd symmetry, but additionally presents the translation invariance along x = z. The curl-free vector field is invariant along the z-axis and equivariant under mirror reflections across the coordinate planes, L i,j : (x, y, z) → (-1) i x, (-1) j y, z , ∀i, j ∈ {0, 1} .This HHD not only uncovers more knowledge of the Chua circuit's symmetry but can also be used to analysis the energy of the system. According to the generalized Hamiltonian formalism (Eq. 20), the energy function H associated with the Chua circuit satisfies the following PDE:EQUATIONwhich is satisfied by the quadratic form:EQUATION)As presented by [77] , this energy function indicates that the Chua circuit keeps oscillatory when the energy release along the x-axis is enough to balance the energy pumping along the y-axis and z-axis.RMSE and MNLL of state derivatives To evaluate the models' performance in terms of learning ẋ = f (x), we first compute the root mean squared error (RMSE) between the predicted time derivatives xi and the ground truth ẋi , i.e.EQUATION)where ∥•∥ 2 is the Euclidean norm, and m is the number of test data. The lower the RMSE, the better. The test set {(x i , ẋi )} m i=1 was generated by sampling a grid of points with the resolution of 10. In addition to the RMSE for evaluating the regression results, we further use the mean negative log likelihood (MNLL) to evaluate the prediction uncertainty provided by the GP models, and it is calculated byEQUATIONwhere var (x i ) is the prediction variance at a test state x i . The lower the MNLL, the more effectively the forecast uncertainty reflects the prediction error.VPT of state trajectories To obtain a more comprehensive evaluation of the learned ODE models, we considered additional metrics for evaluating the accuracy of predicting state trajectories. Recent studies [48, 30, 70] have suggested that RMSEs of state trajectories over long time horizons can be misleading indicators. So we used the valid prediction time (VPT) to measure the model's ability to do such long-term extrapolation in a phase space. Following the definition in [30, 70] , VPT calculates the first time step t at which the normalized root mean square error (NRMSE) between the predicted state xt and the ground truth x t exceeds a given threshold ϵ,EQUATION)whereNRMSE (x t , x t ) = (x t -x t ) T Σ (x t -x t ) /n 1 2, Σ = diag (1/σ 1 , . . . , 1/σ n ), with σ i denoting the variance of the i-th dimension of the true trajectory. The VPT measures how long the predicted trajectory remains close to the true trajectory, so the higher this indicator, the better. In the experiments, we set ϵ = 0.01. To alleviate the dependency on the initial condition, we reported the VPT averaged over trajectories simulated from 50 randomly sampled initial conditions. And a trajectory from an initial condition was solved by the Dormand-Prince method (dopri5) [15] implemented in torchdiffeq4 , integrating forward in time at a frequency of 25 Hz for 15 seconds, with the relative and absolute tolerances of 10 -6 .The experiments were performed on a single Nvidia GeForce GTX 3050 Ti GPU, and all of the models were implemented with PyTorch [51] . The GP-based models (Ind-GP, GIM-GP, Div-GP, HHD-GP and SPHHD-GP) were trained by maximizing the log of their marginal likelihood:EQUATIONwhere |•| computes the determinant of the covariance matrix K + Σ. Training a GP model refers to optimizing the kernel parameters in K and the noise variances in Σ, and these hyperparameters were initialized randomly in our experiments. The GP-based models were trained by the ADAM optimizer [34] , with a learning rate 0.01 for 3000 gradient steps.The Ind-GP models each dimension of a n-dimensional dynamical system independently, so its matrix-valued kernel is constructed byEQUATIONwhere κ i , i = 1, . . . , n are independent scalar kernels. And a standard choice for κ i is the squared exponential (SE) kernel,EQUATIONwhich has two parameters: σ determines the variation of function values from their mean, and l controls the length scale on which the function varies. Realizations of GPs with SE kernels are dense in the set of smooth functions C ∞ (R n , R) (cf. Prop.1 in [37] ). For a fair comparison, we also used the SE kernel (Eq. 55) to construct kernels for the other GP-based models.The GIM-GP produces predictions with the desired symmetry. Each of the systems in our experiments exhibits the odd symmetry as a whole. So their symmetry groups can be given by G odd = {I n , -I n }, where I n is the n-dimensional identity matrix, and the group elements are linear representations of the group actions, i.e., L g (x) = gx, ∀g ∈ G. Therefore, according to Eq. 17, the GIM-kernel for G odd was constructed byEQUATIONwhere the SE kernel (Eq. 55) was used as the basis kernel κ in Eq. 17.GPs with a div-free kernel can be used to approximate conservative dynamics. According to Eq. 12, a two-dimensional div-free kernel was constructed byEQUATIONwhere κ u in Eq. 12 is denoted by κ H instead, indicating that it is the kernel for the Hamiltonian functions of the damped mass spring and the damped pendulum, and the SE kernel (Eq. 55) was used for κ H . The partial derivatives in Eq. 57 were calculated by automatic differentiation in PyTorch, so we did not need to derive its analytic expression.Then we constructed the div-free kernel for the Chua circuit system. According to Eq. 11, a threedimensional div-free vector field is given by f div (x) = Ψu (x), where Ψ = [ψ 12 , ψ 13 , ψ 23 ] ∈ R 3×3with its components given byEQUATION)EQUATION)EQUATION)Then by assuming u = [u 12 , u 13 , u 23 ] T ∼ GP (0, κ u = κ ind ), we constructed the div-free kernel for the Chua circuit according to Eq. 12:EQUATIONEQUATIONwhere κ 1 , κ 2 and κ 3 are all independent SE kernels (Eq. 55).The HHD-GP consists of two independent GPs added together, modeling curl-free and div-free dynamics respectively. GPs with div-free kernels have been constructed earlier, so here we build curl-free kernels according to Eq. 8. For the damped mass spring and the damped pendulum, their two-dimensional curl-free kernel was constructed byEQUATIONSimilarly, the three-dimensional curl-free kernel for the Chua circuit was constructed byEQUATIONwhere κ V is a SE kernel.According to Theorem 5.2, a G-equivariant curl-free vector field is constructed by constraining its scalar potential V to be G-invariant. Therefore, a G-equivariant curl-free kernel is obtained by constructing its potential kernel κ V according to Eq. 16.The curl-free vector field of the damped mass spring and the damped pendulum has two types of symmetry:(i) translation along the q-axis, i.e., L g (x) = x + (g, 0), for all g ∈ R;(ii) rotation by an angle of π radians, i.e., L g (x) = gx, for all g ∈ {I 2 , -I 2 }.These symmetry groups were enforced to a SE kernel sequentially according to Eq. 16, where |G| can be ignored. Specifically, the translation invariance was enforced using the Gaussian integral formula:EQUATIONEQUATIONthen based on which the π rotation invariance was enforced byEQUATIONThe curl-free vector field of the Chua circuit (Eq. 46) has similar two types of symmetry:(i) translation along the z-axis, i.e., L g (x) = x + (0, 0, g), for all g ∈ R;(ii) mirror reflections across the coordinate planes, i.e., L i,j : (x, y, z) → (-1) i x, (-1) j y, z , for all i, j ∈ {0, 1}.Therefore, κ V for the Chua circuit was constructed byEQUATIONG.5.2 Symmetry-preserving div-free kernels Theorem 5.3 shows that a G-equivariant div-free vector field is constructed from a vector potential v with the same equivariance. For the div-free vector field of the damped mass-spring system with SO(2)-equivariance, its κ v was constructed by the GIM-kernel (Eq. 17):EQUATIONThen this κ v (Eq. 70) is substituted into Eq. 18 to construct κ H for its div-free kernel κ div (Eq. 57):EQUATIONEQUATIONwhich admits no closed-form solution, so we used a numerical approximation of the integral (Eq. 72) by sampling discrete rotations of θ = π 4 n | n = 0, . . . , 7 , whose rotation matrices form a finite group. And by setting θ ∈ {0, π}, Eq. 72 was used to construct κ H for the div-free vector field of the damped pendulum, which has the equivariance under a π rotation (odd symmetry).The div-free vector field of the Chui circuit (Eq. 46) has the odd symmetry and the translation symmetry along x = z. κ v respecting the translation symmetry was constructed byEQUATIONEQUATIONwhere p = 1 √ 2 (x + z) , y T . And we further enforced the odd symmetry byEQUATIONThen, κ v = κ • I 3 was substituted into Eq. 18 to construct κ u byEQUATIONFinally, this κ u (Eq. 79) was substituted into Eq. 12 to construct the symmetry-preserving div-free kernel for the Chua circuit.We used the released code5 of D-HNN and ran their training routine for our systems.G.7 Implementation details of predicting energy G.7.1 Predicting energy of the mass-spring system and the pendulum When HHD-GP and SPHHD-GP are used to learn the damped mass-spring system and the damped pendulum, the potential kernel κ H for constructing their div-free kernels (Eq. 57) can be interpreted as placing a GP prior on the Hamiltonian function. Therefore, a joint GP describing both the Hamiltonian function H (x) and the dynamical system f (x) is given byEQUATIONwhereEQUATIONEQUATIONAfter training the model f (x) ∼ GP (0, κ hhd (x, x ′ )) on noisy observations Y = [y 1 , . . . , y m ] T 6 at states X = [x 1 , . . . , x m ] T , we are interested in predicting the value of Hamiltonian function H (x * ) at a new test state x * . Since these data determine H (•) only up to an additive constant, we assume that we have an anchor point H (x 0 ) which can be chosen arbitrarily. Then, according to the GP prior (Eq. 80), H (x * ) andY H = [H (x 0 ) , Y]T are jointly distributed asEQUATIONwhereEQUATIONEQUATIONThen, we obtain the posterior distributionEQUATIONwhere the mean function is used for energy prediction.[60] stated that the div-free component of a chaotic system is always orthogonal to the gradient of its energy function, i.e., ∇H T f div (x) = 0, for all x ∈ R n . The energy function of the Chua circuit is in a quadratic form:EQUATION)So we parameterize the energy function of the Chua circuit by Ĥ (x) = 1 2 a 1 x 2 + a 2 y 2 + a 3 z 2 , where a = [a 1 , a 2 , a 3 ] are parameters. Then by learning a divergence-free vector field fdiv (•) through our model, we can estimate the parameters a by minimizing where Q = q T q with q = fdiv (x 1 ) , . . . , fdiv (x m ) T ∈ R m×3 . And in the experiments, we solved the QP using the solver provided by CVXOPT [2] .The two constituent kernels of our model define the space of divergence-free vector fields (f div ∈ F div ) and the space of curl-free vector fields (f curl ∈ F curl ), respectively. These two spaces overlap partially due to the presence of harmonic vector fields (f harm ∈ F div ∩ F curl ). To eliminate this overlap and make HHD unique, we propose to impose symmetry constraints on the two spaces separately, with the corresponding symmetry groups defined as G div and G curl . Therefore, the uniqueness property of HHD depends on the space of harmonic vector fields that respects the union of two symmetry groups, i.e.F harm = f harm | f harm • L g = J Lg f harm , ∀x ∈ R n , g ∈ G div ∪ G curl .(88) If F harm = ∅, the HHD is unique, i.e. our model is identifiable. A harmonic vector field can always be represented by the gradient field of a harmonic function h, which is a scalar function satisfying ∇ • ∇h = 0 (Laplace's equation). For the three dynamical systems used in our experiments, it can be easily verified that the presence of harmonic components can be eliminated through forced symmetries. The symmetry group union of a damped mass-spring system consists of a rotation group G div = SO (2) and a translation group G curl = {(g, 0) | g ∈ R}. A harmonic vector field f harm = ∇h that satisfies this translation symmetry implies that its harmonic function h (q, p) is independent of the variable q, so the harmonic vector field is given by f harm = (0, ∂ p h) and ∂ p h should be a constant to satisfy the Laplace's equation, but the harmonic vector field in the form of constant clearly contradicts rotation symmetry G div . Therefore, there is no harmonic vector field that respects both the symmetry groups G div and G curl . Similar conclusions can be drawn for the damped pendulum and the Chua circuit. The Laplace's equation and translation symmetry imply that harmonic vector fields can only exist in the form of constant vector fields. However, constant vector fields obviously contradict odd symmetry or mirror symmetry.Compared to the diagonal kernel that independently models each dimension of a dynamical system [13, 31, 68, 69] , our kernel can capture correlations among output dimensions because its off-diagonal elements are non-zero. However, this advantage comes together with a low computational efficiency. The exact inference of HHD-GP over an n-dimensional system with m data points has a cost of O m 3 n 3 in time, due to the inversion of its full covariance matrix K * ∈ R mn×mn . Although in this work we do not focus on the computational challenge in large-scale applications, one can plug our kernel in a variety of sparse GP frameworks [62, 66] , reducing the computational complexity to O M 2 mn , where M (M ≪ m) is the number of inducing point to approximate the covariance matrix.To further compare the predictive performance of the models, we evaluated them by increasing the amount of training data. The results are shown in Fig. 5 . As expected, the performance of the models improves as the amount of training data increases, but we can observe that our model (SPHHD-GP) performs best at every data amount relative to the baseline models, verifying the robustness of our model against data sparsity. J.2 Trajectory prediction of the mass-spring system and the pendulum Fig. 6 and Fig. 7 present the trajectory predictions for the damped mass-spring system and the damped pendulum, respectively.Fig. 8 presents predicted decompositions of HHD-GP and SPHHD-GP. It can be observed that although they capture highly similar system dynamics (the first column in Fig. 8 (a)), they learn completely different decompositions (the second and third columns in Fig. 8(a) ). Compared with the ground truth in Fig. 1 , SPHHD-GP learns the physically correct decomposition, so it can accurately predict the system energy. From Fig. 8 (b) we can find that the predictions of HHD-GP have large variance, meaning that the model is less certain in isolating individual effects from other terms.To demonstrate the interpretability of the learned curl-free features, we adapted the trained model to predict dynamics with different friction coefficients. As detailed in Appendix A.1, the curl-free dynamics in the mass-spring system and the pendulum is caused by friction forces in the system. Utilizing this interpretation and the additive structure of the HHD-based models, we can generalize the models to perform inference over dynamics with different friction coefficients. We performed this experiment by first training HHD-GP, SPHHD-GP and D-HNN with data when ρ = 0.1 (ρ is the friction coefficients), then generalizing the trained models to dynamics when ρ = 0.05 and ρ = 0.5 by multiplying the learned curl-free component with constants 0.5 and 5.0, respectively. The results are shown in Fig. 9 , where we can observe that the three models accurately predict the trajectories of the systems with a friction coefficient of 0.1 (first column in Fig. 9 ), but HHD-GP and D-HNN are difficult to generalize to cases with friction coefficients of 0.05 and 0.5 (second and third columns in Fig. 9 ). In contrast, SPHHD-GP effectively captures the dynamics of different friction conditions, even though it has not been trained for these specific friction coefficients. This adaptability is valuable because it allows the model to be applied to real-world scenarios where friction coefficients may change.To evaluate the robustness of the models, we added the experiments of increasing the standard deviation σ = {0.01, 0.05, 0.10, 0.20} of Gaussian noise in the training data. Fig. 10 presents the results on the pendulum. As expected, the performance of the models degrades as noise increases. However, it can be observed that SPHHD-GP still performs best at all noise levels compared to the baselines, meaning that our model is more robust to data noise. Furthermore, from Fig. 10(d ) we observe that the energy predictions of HHD-GP and D-HNN generate consistently significant errors, insensitive to noise level in training data, which again demonstrate the non-identifiability suffered by these models. In contrast, the energy errors of SPHHD-GP exhibit an increasing trend with increasing Generally, our model can have a significant impact on a wide range of applications, such as physical simulation and robotic control. However, it is imperative that we meticulously evaluate the performance of prediction and uncertainty estimation when implementing it in societal contexts.NeurIPS Paper Checklist Justification: High computational complexity is the main limitation of our proposed model, which has been discussed in Appendix I. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.m is the number of entries above the diagonal. Each off-diagonal element of the matrix corresponds to a scalar function, with elements below the main diagonal as the negatives of those above.For G ⊆ E (n), it holds that ∥Lg (x) -Lg (x ′ )∥ = ∥x -x ′ ∥, for all x, x ′ ∈ R n , and g ∈ G. Therefore, κ (x, x ′ ) = κ (Lg (x) , Lg (x ′ )) is satisfied if κ is an isotropic kernel, i.e., κ (x, x ′ ) = κ (∥x -x ′ ∥), common examples of which are the squared exponential kernel and the Matérn class of kernels (cf. chap.4 in[53]).https://github.com/rtqichen/torchdiffeqhttps://github.com/DrewSosa/dissipative_hnnsyi = f (xi) + ϵ, ϵ i.i.d ∼ N (0, σI) 


Decision: Reviewers agree that this manuscript on gaussian processes kernels for the Helmholtz-Hodge Decomposition of dynamical systems is clearly written and novel. This is an exciting contribution with many potential uses.

Reviews: review: paper_summary: paper_summary: This paper derives a G-equivariant versions of the both curl (through the Haar Integration kernel) and divergence free (through GIM kernel) kernels. This is then used to define a prior over the Helmholtz decomposition and is identifiable wrt to relevant functions in the Euclidean group (translation, rotation, reflection etc). This is demonstrated on a range of examples.
strength_weakness: strength_weakness: 1.	The justification for how and why the curl and div free kernels should be invariant is not given until the experiments/the appendix. It would be nice if this was discussed in section 5.
2.	The fact that the actual symmetries required is problem specific is not really clear in the main paper. In the appendix it is stated, but it would be good if this was made more explicit in the main paper.
questions: questions: 1.	Is the proposed symmetry-preserving Helmholtz decomposition as expressive as non-symmetry preserving one? Ie any (sufficient smooth) field can be decomposed into curl and divergence free part. Is the same for the symmetric-preserving decomposition?
2.	On line 328 it is stated that mirror symmetry is incorporated into the div-free kernel. Where is this kernel given?
3.	Are there general symmetries that should be held or is this always going to be problem/application specific? 
4.	When the kernels do not admit a closed form how sensitive are they to approximations? And does this significantly affect the computational cost?
5.	Is there any impact of using sparse GPs? Ie are the sparse models still guaranteed to be G-invariant? 
6.	Relevant GP citations:
    o	Equivariant Learning of Stochastic Fields: Gaussian Processes and Steerable Conditional Neural Processes, Holderrieth et al
    o	Vector-valued Gaussian Processes on Riemannian Manifolds via Gauge Equivariant Projected Kernels – Hutchinson et al
limitations: limitations: yes
paper_summary: paper_summary: The paper formulates a vector-valued GP model to infer unknown vector fields. The authors discuss how to impose symmetry-based constraints on the GP models to ensure physically meaningful decompositions. The paper provides theoretical proofs to support the construction of curl-free and divergence-free vector fields that preserve desired symmetries. The model's effectiveness is validated through experiments on various dynamical systems, including dissipative Hamiltonian dynamics and chaotic systems.
strength_weakness: strength_weakness: - There are missing citations for relevant existing methods.
- There is no discussion on computational complexity.
questions: questions: - In [R1], a Gaussian Process for dissipative Hamiltonian systems has been proposed. Can HHD-GP be considered a generalization of [R1]? What are the essential differences between SSGP and HHD-GP? While SPHHD-GP is clearly novel, it was unclear whether HHD-GP is new compared to [R1].
- Please add a discussion on computational complexity. Specifically, numerical integration is used to consider symmetry constraints, but could this be a computational bottleneck?
- In Figure 3, the error of HHD-GP appears to increase with the number of samples, which seems unnatural. Shouldn't this be analyzed in more detail?

[R1] Yusuke Tanaka, Tomoharu Iwata, Naonori Ueda, Symplectic Spectrum Gaussian Processes: Learning Hamiltonians from Noisy and Sparse Data, NeurIPS, 2022.
limitations: limitations: The differences between HHD-GP and existing methods (such as [R1]) should be thoroughly discussed.
paper_summary: paper_summary: The authors tackle the problem of modelling dynamical systems in scenarios where it may not be possible to straightforward to determine the exact form of the ODEs governing the system and optimise their parameters directly. Whilst physics-informed Bayesian models which learn divergence-free vector fields are capable of extrapolating dynamics effectively in an interpretable manner, they cannot represent certain common real-world behaviours such as dissipation. Proposed in this work is a Gaussian process (GP) which combines the utility of the div-free vector field with a curl-free vector field, an approach which is theoretically motivated by the Helmholtz-Hodge decomposition which states that many vector fields may be decomposed into a sum of a div-free and curl-free field. This approach (termed the SPHHD-GP) allows for a much wider range of dynamical systems to be modelled without strict assumptions on the form of the system, and also permits a symmetry preserving representation which allows for enhanced interpretability of the learned dynamics. Experimental results are presented which show the efficacy of the approach on tasks such as learning Hamiltonian dynamics from noisy data, learning chaotic dynamics and modelling ocean currents.
strength_weakness: strength_weakness: -	From the text it appears the ocean current experiment was performed using a sparse GP implementation, whereas the rest of the experiments were performed using exact inference, is this correct? This is only mentioned very briefly in passing and few details given; I know the details of the sparse GP formulation are probably a given to readers with a background in the area, but I think adding some detail to the appendix on how this was implemented would be useful to many readers.
-	Several of the derivations, proofs and discussions in the appendix are of sufficient relevance that it would be great if they could be included in the main text, however I’m obviously aware that the authors are constrained by the page limit so I wouldn’t expect any alterations based on this. I think the amount of content more just speaks to the fact the work has been carried out in a thorough manner.
questions: questions: -	Please address the point mentioned in the weaknesses section regarding adding some brief details about the sparse implementation into the appendix and make some reference to this in the main text (maybe where the sparse GP framework is mentioned in Section 6.3).
-	Related to the above point, could you speak to how the sparse formulation performs empirically compared to the exact approach on problems for which both are computationally feasible (i.e., sections 6.1 and 6.2)?
limitations: limitations: The authors address the main limitation of the model in Appendix I which is the considerable computational complexity of the model, $\mathcal{O}(m^3 n^3)$, although this approach can also be implemented using sparse GPs instead to reduce this. The authors also briefly mention in Appendix K how the model has a wide range of applications, but from a societal impact perspective it is important that predictions and uncertainty estimates are rigorously evaluated. As I don’t think the work has significant potential for negative impact, I think this is adequate.
paper_summary: paper_summary: This paper proposes a novel dimensionality reduction method.
The method relies on the Helmholz-Hodge decomposition to identify the dynamical system through a decomposition into a curl-free and a divergence-free part.
Furthermore, the method introduces a way to incorporate priors that constrain the identified vector fields based on symmetry constraints.
This method is applied to various toy models (mass-spring system, pendulum, Chua circuit) and to real world data (ocean current fields).
strength_weakness: strength_weakness: The overall applicability of this method is questionable due to not completely addressing the identifiability issue of the decomposition. The issue is only resolved for systems with known symmetries, see also Questions.

The paper mentions that the accurate prediction of energy demonstrates the interpretability of learned div-free features. 
However, this is not well-explained. The paper claims that the curl- and div-free components are physically meaningful, but this should be better justified.

The paper is missing comparison to other existing approaches for modeling physical systems. See for example: 
Guo, Q., Mandal, M. K., \& Li, M. Y. (2005). Efficient Hodge–Helmholtz decomposition of motion fields. Pattern Recognition Letters, 26(4), 493-501.

Compare to other dynamical systems reconstruction methods:
- Hess, F., Monfared, Z., Brenner, M., \& Durstewitz, D. (2023). Generalized teacher forcing for learning chaotic dynamics. arXiv preprint arXiv:2306.04406.

-I. Mac\^edo and R. Castro, “Learning Divergence-Free and Curl-Free Vector Fields with Matrix-Valued Kernels,” technical report, Instituto de Matema ́tica Pura e Aplicada, Rio de Janeiro, Brazil, 2010.

Finally, there is no code provided to show the method's working and not all hyperparameters are described (see also Questions).
questions: questions: What is the relevance of $\mathcal{G}$-equivariance and uniqueness/identifiability?
Appendix H: Uniqueness and Symmetries
- How do you know what priors to choose? Isn't an agnostic method better?
- How do you know which symmetries to enforce?

Could it be more fruitful to restrict yourself to a compact manifold with boundary so that the decomposition is unique?



# Implementation details
- How were the kernel parameters and noise variances initialized? Does training depend on initialization?
(Here you mention that these are hyperparameters, but typically one refers to non trained parameters as hyperparamters.)
-What were the ADAM optimizer paramters used for training?


# Choice of evaluation metrics 
-Why did you mostly focus on performance of the vector field? If the point is learning the dynamics, metrics that look at the dynamics directly would be preferable.
Why not prediction of the future state/trajectory? Why prediction only measured through VPT?
limitations: limitations: The authors should comment on limitations and societal impact.

Comments: title: Rebuttal response
comment: Thank you for your reply. My concerns are now largely resolved, and the differences between the proposed model and SSGP are much clearer. I will raise my score from 5 to 6.
Thank you for your thorough response. I think including the discussions above in the paper will improve an already well-written paper.  I have no more questions and maintain my original score of acceptance.
Thank you for your rebuttal. However, I still believe that the paper would benefit from comparisons to relevant baselines (besides GPs and HNNs). Specifically, including a comparison to other dynamical systems reconstruction methods would provide valuable context and help to clarify the advantages of your approach. Such comparisons could significantly strengthen the paper by demonstrating its effectiveness relative to existing techniques. 

I will keep my score as is.
Thank you very much to the authors for taking the time to clarify those points, I appreciate it's primarily details that will be apparent to people closely familiar with GPs, but I do think expanding on them will broaden the accessibility of the paper to those coming more from a physics background.

I keep my score unchanged and recommend acceptance.