HHD-GP: Incorporating Helmholtz-Hodge Decomposition into Gaussian Processes for Learning Dynamical Systems
Authors: Hao Xu
Jia Pan
Venue: NeurIPS 2024 poster (2024)
Keywords: Gaussian process
Helmholtz-Hodge decomposition
dynamical system
Primary area: machine_learning_for_physical_sciences

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.

HHD-GP: Incorporating Helmholtz-Hodge Decomposition into Gaussian Processes for Learning Dynamical Systems
Hao Xu
2 Jia Pan 1.2,3+
'Te University of Hong Kong, Hong Kong, China
Centre for Transformative Garment Production, Hong Kong, China
LimX Dynamics
xuhaovic@connect hku. hk, jpan@cs hku.hk
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-agnosia that lack interpretabilltyy robustness, and data efficiency. This paper mitigates this concern by incorporating the Helmholtz-He Hodge decomposition into Gaussian process model, leading to a versatile framework that simultaneously learns the curl-free and divergence-free components of dynamical system. Learning predictive model in this form facilitates the exploitation of symmetry priors. In addition 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.
design, leading learning methods
Introduction
dynamical system describes how the state of system evolves over time (64. Data-driven modeling of dynamical systems has become a fundamental task in many modern science and engineering applications, such physical emulation (21) and robotics control [33. Mathematically, dynamical system is often a set of first-order ordinary differential equations (ODEs) or, equivalently, a smooth vector field on 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 and non- parametric Bayesian models 663 27 26 for the black-box approximation of vector fields. Despite the rich expressive power, the when applied physical systems. To address this issue, 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 learning Hamiltonian systems 23 67 54] and incompressible fluid dynamics [73 33),
design of these models hinders their performance
(div-free) vector fields, as 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
ML models are constructed to learn divergence-Fif
"Corresponding author.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

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, explore supplementing the div- free vector field with a curl-free vector field. This is inspired by the Helmholiz-Hodge decomposition (HHD)
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-frre vector field. HHD is widely used in the study of Navier-Stokes
3
equations 45 20 91, but in this work we explore its connections with more general dynamical systems. For example, as shown in Fig. [ HHD can be used to characterize the dynamics of dissipative Hamiltonian system, where the div-free component represents its conservative dynamics and the curl-free component describes the friction friction-induced induced dissipation. In addition to dissipative Hamiltonian systems, HHD 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
In this work, we construct Gaussian pro- cess (GP) priors for div-free and curl-free dynamics separately, resulting in an addi- tive GP model in the form of HHD for learning vector fields of dynamical systems. Our resulting HHD-GB allows one to lever- age more prior knowledge than modeling the system as a whole. In particular, we in- vestigate its potential in exploiting priors of symmetries, motivated by the observation that the div-free and curl-free components of dynamical system usually exhibit more symmetry properties than the system itself.
(a) dynamics
(b) div-free part
(c) curl-free part
Figure I: The vector field in phase space of mass- spring system and its HHD.
For example, the damped mass-spring system in Fig. [(a) exhibits odd symmetry, but its div-free (Fig. [(b) and curl-free (Fig. [[c] components additionally present rotation and translation symme- try, respectively. Therefore. 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-meant menniigful............ In particular, the learned div-free features are closely related to the energy of dynamical systems.
prrseetty-reserviing
The main contributions this work are summarized as:
We introduce a GP prior (called HHD-GP) for learning the Helmholtz-Hodge decomposition of general dynamical system.
We construct a symmetry-preserring 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.
2
Related work
Learning with div/curl-free constraints Div-free vector fields are 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 universal representation of div-free vector fields. Based on the same representation, we constructed div-free kernels for GPs from matrix- caatrixvaaledd 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 77 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 ours when the dimension of HHD-GP is two, but we further developed
extension of HHD-GP to solve its non-identifiability problem. Another similar work Dissipative Hamiltonian neural network (D-HNN) [22), which compensated HNN (23] with 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.
pymmeetr--reseervee
Learning with symmetry Symmetries are another important aspect of priors that can be incorpo- rated into machine learning models. Motivated by the success of the translation-invariae NNs (38), network architectures with symmetries to more general transformations have been proposed, such as steerable CNNs [74] [O] and graph NNs 46 6 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. [4] 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 impose Euclidean symmetries div/curl-free vector fields, which to the best of our knowledge has not been explored by the machine learning community.
3 Background
3.1 HHD and Problem Setup
We consider an autonomous system governed by the following ODEs:
* (t) = dx (t) f (x (1)) = fourt (x (1)) + faio (x (1)),
di
which defines vector field by assigning vector (x) E R" to every state x E R". We assume that the vector field f E 1.2 (RR", R²)) is smooth, and any such vector field can be decomposed into the sum of a curl-free vector field fourt R. - Ih" (V A feurt 0, Vx E IR") and divergence-free vector field fiie= R - R. (V faiv 0, vx E RR),, according to the Helmholtz-Hodge decomposition (HHD) [44 88. In this work, we are interested in learning f, fownt and faiv simultancously from collection of noisy observations denoted by D = {(x,, y.) ) with a noisy observation y; at a state x; given by
y: f (x.) +,,
i.d ~ (0,
(2)
where an additive noise e E R2 follows a zero-mean Gaussian distribution defined by a covariance matrix Q = diag (o?...., o E Rxx modeling noise variance in each output dimension.
3.2 Vector-valued GP Model
We are interested in using Gaussian processes (GPs) to infer unknown vector fields. A GP is a stochastic process commonly used 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" + R", we assume
valued GP prior as f (x) GP (0, x xx)x where the mean of the function values is set to zero, and the covariance is captured by a matrix-valued kernel A : RR x IR" + IR"x", whose (i, j-th entry expresses the correlation between the i-th dimension of (x) and the j-th dimension of f (x'). In the GP framework, the kernel controls the properties of possible functions
under GP, leading various efforts of problem-spectee design of kernels [76 [6 18] And to be valid covariance function, the kernel should be symmetric and positive semidefinite 788
GPs provide Bayesian nonpparametric approach for solving regression tasks. According to the GP prior, function values at inputs X = [xi.... x,,]' are jointly distributed as M (f (X); O. K), where K [K xx) E Raxn m kij=1 is block bookk-pattiiioned covariance matrix. Then the marginal
likelihood of the noisy observations Y [y1...., yml T given by Eq. 2 be calculated by
p (Y (Y / (x X) = N (Y O.K.),
(3)
3

---

where K. = - & Im is a diagonal matrix whose elements are the variance of the observation noise. Training a GP model refers maximizing the log of 3 optimize the kernel parameters and the noise variance. Then by conditioning on these observations using Bayes" rule, the predictive posterior for new state x, is still Gaussian with its mean ja and variance v given by
# (x.) = k. Y, v (x.) = K (Xx, x.) - k k"K,"' k,
where k = [<: (x1, (x,,))) ..... (Xm, x.)] E R= The derived mean function is used for regres- sion 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.
4 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 fourt and fase, respectively, fourt ~ GP (0, Keur!), fdin -GP (0, N. 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- freee ones. And the additive kernel K = Kicurl + Ediio inherits the symmetric and sosidive-eemidefinite properties of Kear! and Kidivs 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. curl-free part Mourt and a divergence-free free part Mdivs and we have
4 can be split into a
μ Peart + Pdiv kurrt 'Y + K.
(5)
where K. = Keurt + Kaiv + Σ. It can be seen that the effects of feurt and fate can be treated as observation noises for each other, so their prediction variances at x. are obtained by
Veurt (x.) = Keart (x,,
kT
K K_ ' Vdieo (x.) = Kddiv x.) J., k (6)
Consequently, observations of f (x) given by Eq. 2 can be used make predictions for its hidden components four! and fdiu- 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 Keurt and Kaio from the representations of fours and lai respectively, exploiting closure of GPs under linear transformation. Let R = R ...,...., 8z,.] be the polynomial ring in the partial derivatives, and L E Raxb be a matrix of differential operators acting on functions g R R distributed as GP (μ (x). K x')). Then, the transformation of under L is again distributed as a GP with
Lg GP (x), Lxx (x, x') CT. R R + Baxa),
where Lx and Cx denote the operation of L on the first and second argument of (x, x'), respec- tively [29 [ 36 77 To make Lx (x. x') LT, valid covariance function, its underlying kernel (x, x') must be twice differentiable in R". See Appendix B for details on linear operations on GPs.
4.1 Curl-free Kernel
The gradient operator (V := [a......,, Io,,'' E Rmx1) defines a surjective mapping from the space of smooth scalar fields to the space of curl-free vector fields [14), so fourt can be represented by fourt = VV, where E ( (IK", .R) is called the scalar potential of fours. Since the gradient operation defines a linear transformation, if GP with scalar kernel Kv is assumed on V, the distribution of feurt is again a GP. According to Eq. 7 the curl-free GP over feurt is given by
Fourt GP (0, Keart V V (x,x') V
(8)
where Keart RR" x R + Rexn is matrix-vall valued kernel constructed by the Hessian of the scalar kernel kv, consisting of all second-order partial derivatives x and x', with the entry the i-th row and j-th column given by
[Kewrl (x, x ) = Cov [aV (x) /0z, OV (x) /0a]] = Paky (x, x') /Or,0r'.
By this construction, if kv induces GP with realizations dense in C (R" R), the set of realizations of GP (0, Neurt) dense in the space of curl-free vector fields, because a surjective mapping maps dense sets to dense sets.

---

4.2 Divergence-free Kernel
div-free vector field can be constructed from Akem-eymmettic matrix field [5 32 57 let A IR" + R"x* be skew-symmetri mutrix-valued function, then a div-free vector field fave can be represented by taking row-wise divergence of A, i.e.,
Specifically,
faie
I AI...... V A.
where A, : Rn + R is the i-th row of A. The skew-symmetric matrix field A of size n x a can be compactly represented by m (n 1) / scalar functions Uij E C (R",R).
0
U12
M12
0
U2n
π-
n
A
Pijliit
(10)
1=++1
U2n
where Biy E "n*=" a matrix with its ji, -h entry equal to I, (j,))tth entry equal to -1, and all
other entries equal to O. Then, the div-free vector field given by Eq. can be reformulated as linear transformation:
fdiv (x)
Σ 1=1 j=
Σ
Viguiy (x) Yu (x)
(1I)
where bi 2 E R is column vector obtained by linearly transforming the gradient operator. The column vectors Viy are aggregated in a matrix T E Raxm. and the m corresponding scalar functions Viy are collected in vector-valued function u R + RT. V [-] is matrix of linear differential operators, so to use a GP to model fdie, we can proceed by assuming a GP prior over
GP (0, Kig IR= x IR" + Rmxm) then based on the closure of GPs under linear transformation the GP prior over faiv (Eq. is constructed by
(Eq.
fdiv
GP (0, x (x, x') TT)
(12)
where Ku is scalar-valued kernel for two dimensional systems (n 2, m 1), and is matrix- valued kernel for n > 2. Notice that this div-free kernel can be treated as generalization of the div-free kernel derived from scalar kernel [50] 75), which can be recovered by setting Ku R- I. Therefore, it is more expressive and flexible. Theoretically, the GP model given by Eq. [[2 can be used to approximate arbitrary div-free vector fields, because the representation given by Eq. is shown to be maximally expressive (i.e. universal) (7)
4.3 Identifiability and Constraints
With the curl-free and div-free kernels, objective is to learn physically interpretable representations of a dynamical system based on the HHD-GP model. However, the HHD always not unique due to the existence of harmonic components fharm (vector fields satisfying both V A fharm 0 and fharm 0. e.g., constant vector fields). For the HHD of a dynamical system with the true functional decomposition fount and five
f (fourd + fharm) + (fāro fharm)
(13)
is a valid HHD for arbitrary fearm. 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 fourt. fato are physically meaningful.
To mitigate the identifiability problem in additive regression models, an effective method is to impose 47 42). The imposed constraints affect the decomposition results of the additive models. Therefore, to ensure that HHD--P 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, the next section, we present how to impose symmetry-baseed constraints on the curl-free and div-free GP models.
constraints their component models [16] 17
2m is the number of entries above the diagonal. Each off oif-diagnnal element of the matrix corresponds to a scalar function, with elements below the main diagonal the negatives of those above.
5

---

Symmetry Constraints
5.1 Equivariance and Invariance
Symmetry a fundamental geometric property prevalent dynamical systems natural and usually described by the concept of equivariance and invariance:
Definition 5.1 (Equivariance and Invariance). Let g be group acting R" through smooth map L:gx R" + R". The dynamical system f R + R" is said to be G-equivariant
(f o Lg) (x) - Jr, (x) f (x) Vx E 9 EG
(14)
where Lg (x) == L Xg.x,, and JL, denotes the Jacobian matrix Lg. Then, 9 is termed the symmetry group of the dynamical system. In particular, if J r., is the identity matrix (i.e., f Lg f. V. E 9), 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 's trajectory commutes with the action map. For vector fields on R" the symmetry group is commonly subgroup of the Euclidean group F: (n), which comprises all intuitive geometric transformations in Rn (see Appendix C for brief introduction). The symmetry constraints refer to that 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.
5.2 Symmetry--pesserving Curl-free GP
The curl-free GP (Eq. S constructed by transforming another GP over a potential function, implying that we can impose constraints of symmetry on the curl-free GP by designing 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 be Euclidean group or its subgroup, and let V .R - R. be G-invariant scalar function. Then, the curl-free vector field feurt : TR" - TR" defined by fourt (x) = VV (x) is G-equivariant.
See Appendix D.I for the proof. Theorem 5.2 shows that G-invariant scalar potential V can yield a G-equivariant gradient field, indicating that if any realization V of GP (0, is constrained to be G-invariant, then its pushforward GP over VV ~ GP (o, V V can induce the space of G-equivariant curl-free vector fields.
is obvious that G-invariant scalar potential V can be constructed by integrating some non-invariant function h _ + R over the symmetry group: V V= SG (h o Lg) dg. where the measure dg is called Haar measure, which exists for locally compact topological groups and finite groups. Therefore, by assuming that h distributed as h GP (0, Kr.), we can construct the GP prior over the G-invariant scalar potential as V ~ GP kv), with its kernel Av given by
[/ h (L, (x)) dg, | (L, (L, (x)) dg = 16/ Kh (Lg (x), (x')) dgdg'.
kv =
Cor Cov
(15)
This kernel is called the Haar- interiintegration 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 Ky is invariant to any G in the sense that (x,x') k (Z.g (x), I.g (x' [x') complexity reduction of Eq. by one square-root can be performed by
KV = [ nn (x, (1g-t)))))))) dgdg' = |9l 101/ Kn (x, T.g) dg,
(16)
where [g] = [g dg, and denotes the cardinality of g when the group is finite.
Foor 9 C E (n). holds that ||4 (x) - Lg (x')l = ||x - x'Il. for all x, x' E R", and E G. Therefore. (x,x') = N (L, (x) Ly (x')) is satisfied if x is an isotropic kernel, i.e., x (x, x') x (]]x - x [l), common examples of which are the squared exponential kernel and the Matém class of kernels (cf. chap.4 in (53).

---

5.3 Symmetry-preserve Divergence-free GP
To incorporate the equivariance condition (Eq. [4] into realizations of the div-free GP (Eq. 12 we construct the skew-symmetric matrix field A from vector-valued function. Specifically, given smooth vector field v E C (IR", ,IR",, A is constructed by A = Jv Jy. where Jv denotes the Jacobian of v with its (i. --th entry given by /v,/Ozz. Then the component function Wy in Eq. I0 is given by uij = 00, - 0v,/āx; /dz By this construction, the symmetry of the div-free vector field fdiv is governed by the symmetry of its vector potential In particular, G-equivariant v can produce a G-equivariant fdio, and it is formalized in the following theorem:
Theorem 5.3. Let 9 be a Euclidean group or its subgroup, and let v : IR" - R2" be a G-equivariant vector field. Then the divergence-feee vector field fase- = (V Ai..... V An] is G-equivariant, where A; denotes the i-th row of the skew-symmetric matrix-valued function A Jv J
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,Ki)) and to constrain v to be g-equivariant, we build its kernel Ni _2xx in the form of the Group Integration Matrix kernel (GIM-kerne)) [56 55], which is constructed by:
Ry (x,x') - f 1xxxx (x, by (x')) J dg,
(17)
where x. is some arbitrary scalar-valued kernel satisfying (x,x') K (Lep (x), Lg (x') for all @ E G. The GIM-kernel spans Reproducing Kernel Hilbert Space (RKHS) of functions with the desired equivariance [55). So we can then use Ky (Eq. 17 to construct the GP prior over u in Eq.
where the covariance between components May and ukg given by
Vig = au Ouy
Ouk aug
[Kulig,tg = Cov
4ky
Or,
82
8²
82
8² dr,ūsr, [ivljxx
[Nvle, + Or,
[
[Nwlin
(18)
Or,
Our,Oorrk
Finally, this matrix-valued kernel Nu 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.
6.1 Learning Dissipative Hamiltonian Dynamics
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,)) (xx 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 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--peservinng preserving extension, SPHHD-GR 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 sealar GP. See Appendix G for the implementation details of these models.
7

---

Table I: Comparison of our models to baselines. The results are averaged over 10 independent experiments performed by resampling the training sets and model initial parameters. The RMSE and the VPT are recorded in the scale of x x10-2 and in the form of mean + standard deviation. Bold font indicates best results.
Damped Mass Spring
Damped Pendulum
Model
RMSF L
VPT
MNL...
RMSE
VPT
MNLI
D-HNN
35,58 6.08
0.67
0.13
NAA
183.38 + 27.92
0.39 0.07
NIA
Div-GP
+: 38,11
2638088 + 7917.11
55.21
15.93
.001 0.36
4.15
1.38
-2.25
80.59
38.08
-098
GIM-GP
2.80 + 1.96
5.28
2.09
-2.80
+ 0.45
26.46 +: 13.83
2.32
0: 0.74
-1.15 +: 0.22
HHD-UP (ours)
4.83 + 1.60
2.92 0.90
-2.20
+ 0.29
34.64 + 15.68
1,33
0.42
.0.9 0.26
SPHEID GP (ours)
2,02
1.75
8.52
5.11
-3.16
0.88
13.34
6.89
1.37
-1.56 0:64::00 0.31
The results for 20 training data are shown in Table I and the results for an increasing number of training data are provided in Appendix [.I] 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 HHDDGP better than that of another HHDD-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-GF performs best overall in learning ODEs. Appendix 7.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 's Hamiltonian, as long as the cor- rect form of HHD is learned by our model. To show this, we evaluated our model by another task: we predicted the Hamiltonian H (x,) (en- ergy) along the system's 's trajectory )x,) For HHD-GP and SPHHD-GP the Hamiltonian is predicted from joint GP prior over the Hamil- tonian and its underlying dynamics (see Ap- pendix G.7.I 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 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 [..3 for an visualization of predicted decompositions, which shows that although the models capture similar system dynamics, they can learn completely different decompositions.
Mass-spring System
Pendulum
Smoud Iuin
2.5
0.2
3.0
0.0
25
1.3
1.00
T500
0.0
2.5
5.0
3.5
160
12.5
15.0
Time
Time
Figure 2: Energy prediction along the trajectory initialized (1.0, 0.0) of the mass-spring system (left figure) and (1.5, 0.0) of the pendulum (right figure), respectively.
Consistency is a necessary condition for a learn- ing model to be identifiable and refers the property that its parameter estimates should
Mass spring System
Pendulum
AMSE Energy
RNSSE of Energy
converge the true values as the amount of
S-PAM SHID-6
DelM ShR12-07
data increases [19] 40]. Therefore, to further explore the non-identifiabilit problem, we pro- vide 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 pa- rameters cannot converge as the amount of data increases. In contrast, the energy errors of SPHHD-GP exhibits decreasing trend, which reflects
-0-3
100
250
315
100
200
300
400
Number data
Number of data
Figure 3: RMSE of energy prediction with increas- ing number of training data.
8

---

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-identifiabiiit 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 [144 we added experiments of adapting the learned models 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 the baselines.
6.2 Learning Chaotic dynamics
Learning the correct HHD also has im- portant implications for studying chaotic systems. [60] stated that the div-free component of chaotic system is always orthogonal to the gradient of its energy function, i.e., H faio (x) = 0. This energy function can be used to synchro-
Table 2: Experimental results on the Chua circuit. R.SE and RMSE refer to the errors in derivative and energy predictions, respectively.
Model
RMSK
MNLL
RMSE
Div-GP
122.4 à: 51,4
0,89 0:74
09 4 0.1
N/A
nize two chaotic systems [60], analyze
Inl-GP
12.5 : 12.6
-1.76 4 1.02
59
N/A
their stability 1771, and design energy- modulated controllers [43. This experi-
GIM-GP
10.2 10.6
-1.92 0.67
63 4 2.9
N/A
HHD-Gi (ours)
12.0 11.1
-1.76 +: 0.68
5.1 + 1.4
7.4 3.6
ment aims to learn the Chua circuit [49]. which is a chaotic system with applica- tions in various fields. Its governing equa-
SPHAD-GP (ours)
4.1 1.8
-2.72 + 0.29
13.0 5.6
0.7 1.5
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 [ which shows that
tion and the HHD are detailed in Appendix E.2
SPHHD-GP GP consistently outperforms the other approaches. Appendix [.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.
6.3 Learning ocean current fields
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 '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, used 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 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 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 enforce the parallel boundary condition. Specifically, if vector field has mirror symmetry with respect to 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 the ocean current field. reduce the computational complexity, the kernels were fitted into sparse GP framework [2) which reduced the computational complexity from O (m²) to O (M2m), 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.

---

Preditiion of Ocean Carrents
Oivergence free Compenent (Vortiitt)
Curl free Component (Chvengence)
(a) Rresults the Helmholtz-Gi
Prodiction of Carents
Divergence Compenent: (Vortiiit)
Component (Dlvergence)
(b) Rresults of the SPHHD GP (ours)
Figure 4: Predictions of the ocean current field and its HHD by (a) Helmholtz-GP and (b) SPHHD-GP The black arrows are the predicted currents, and the red arrows denote the observed buoy velocities. First column: Predicted ocean current fields. Second column: Predicted vorticity fields. Third column: Predicted divergence fields.
The predictions of the ocean current field and its HHD shown in Fig. 4 Although the ground truth 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. The oceanographers expect to find rich structure in divergence predictions. However, the Helmboltz-GP fails to recover the divergence field and instead predicts an almost constant curl-free vector field, as shown the third block of Fig. 4(a) This is caused by the non-identifabiiiiy 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.
7 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.
Acknowledgments
This research partially supported by the Innovation and Technology Commission of the HKSAR Government under the InnoHK initiative and GHP/126/2IGD, and GRF 17200924.
10

---

References
[1] Christian Agrell. Gaussian processes with linear operator inequality constraints. Journal of Machine Learning Research, 20:1-36, 2019.
[2] Martin Andersen, Joachim Dahl, and Lieven Vandenberghe. CVXOPT: Convex optimization. Astrophysics Source Code Library, pages ascl-2008, 2020.
[3] George B Arfken and Hans J Weber. Mathematical methods for physicists, 1999.
[4] Vladimir Igorevich Amol'd. Mathematical methods of classical mechanics, volume 60. Springer Science Business Media, 2013.
[5] Cristian Barbarosie. Representation of divergence-free vector fields. Quarterly of applied mathematics, 69(2):309-318 2011.
[6] Renato Berlinghieri, Brian L.. Trippe, David R. Burt, Ryan James Giordano, Kaushik Srinivasan, Tamay Özgökmen, Junfei Xia, and Tamara Broderick. Gaussian processes at the helm(holtz): more fluid model for ocean currents. In Proceedings of the 40th International Conference on Machine Learning, pages 2113-2163, 2023.
[7] Andreas Besginow and Markus Lange- HngeeHHee Constraining Gaussian processes to systems of linear ordinary differential equations. Advances Neural Information Processing Systems, 35:29386-29399, 2022.
[8] Harsh Bhatia, Gregory Norgard, Valerio Pascucci, and Peer-Timo Bremer. The Helmholtz- Hodge decomposition- survey. IEEE Transactions on visualization and computer graphics, 19(8):138-1140) 2012.
[9] Jean-Paul Caltagirone. On a reformulation of Ntvier-Stokes equations based on Helmholtz- Hodge decomposition. Physics of Fluids, 33(6), 2021.
[10] Gabriele Cesa, Leon Lang, and Maurice Weiler. A program to build E(N)-equivariant steerable CNNs. In International Conference on Learning Representations, 2021.
[11] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[12] Alexandre Joel Chorin, Jerrold E Marsden, and Jerrold E Marsden. A mathematical introduction fluid mechanics, volume 3. Springer, 1990.
[13] Marc Peter Deisenroth, Dieter Fox, and Carl Edward Rasmussen. Gaussian processes for data-efficient learning in roboties and control. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2):408-423, 2015.
[14] Manfredo P Do Carmo. Differential forms and applications. Springer Science & Business Media, 1998.
[15] JR Dormand and PJ Prince. A reconsideration of some embedded runge kutta formulae. Journal of computational and applied mathematics, 15(2):203-211, 1986.
[16] Nicolas Durrande, David Ginsbourger, and Olivier Roustant. Additive covariance kernels for high-dimensional Gaussian process modeling. In Annales de la Faculté des sciences de Toulouse: Mathêmatiques, volume 21, pages 481-499, 2012.
[17] Nicolas Durrande, David Ginsbourger, Olivier Roustant, and Laurent Carraro. ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis. Journal of Multivariate Analysis, 11 115:57-67, 05:57-67,20111 2013.
[18] David Duvenaud. Automatic model construction with Gaussian processes. PhD thesis, Univer- sity of Cambridge, 2014.
[19] Arne Gabrielsen. Consistency and identifiability. Journal of Econometrics, 8(2):261-23.. 1978.

---

[20] Abdou Garba and Pierre Haldenwang, A Srumholtz-Hoogee Projection Method Using an Iterative Gauge Computation Solve the 3D Generalized Stokes Problem. SIAM Journal on Scientific Computing, 35(3):A1560-A1 1560-A1583, 2013.
[21] Roman Garnett, Shirley Ho, and Jeff Schneider. Finding galaxies in the shadows of quasars with Gaussian processes. In International Conference on Machine Learning, pages 1025-1033. PMLR, 2015.
[22] Sam Greydanus and Andrew Sosanya. Dissipative Hamiltonian neural networks: Learning dissipative and conservative dynamics separately. arXiv preprint arXiv:2201. 10085, 2022.
[23] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks, Advances neural information processing systems, 32, 2019.
[24] Bernard Haasdonk, A Vossen, and Hans Burkhardt. Invariance in kernel methods by Haar- integration kernels. In Image Analysis: 14th Scandinavian Conference, SCIA 2005, Joensun, Finland, June 12-22. 2005. Proceedings 14, pages 841-851. Springer, 2005.
[25] Marc Härkönen, Markus Hegge-Hngnmm,n, and Bogdan Raitä. Gaussian process priors for systems of linear partial differential equations with constant coefficients. arXiv preprint arXiv:2212. 14319, 2022.
[26] Pashupati Hegde, Çağatay Yıldiz, Harri Lähdesmäki, Samuel Kaski, and Markus Heinonen. Variational multiple shooting for Bayesian ODEs with Gaussian processes. In Uncertainty in Artificial Intelligence, pages 790-799, PMLR, 2022.
[27] Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ODE models with Gaussian processes. International Conference Muchine Learning, pages 1959-1968. PMLR, 2018.
[28] Morris W Hirsch, Stephen Smale, and Robert L Devaney. Differential equations, dynamical systems, and an introduction to chaos. Academic press, 2012.
[29] Carl Jidling, Niklas Wahlström, Adrian Wills, and Thomas B Schön. Linearly constrained Gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
[30] Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis. Sympnets: Intrinsic structure- structure-pesserving symplectic networks for identifying Hamiltonian systems. Neural Networks, 132:166 10211666177 166179,6666611 2020.
[31] Sanket Kamthe and Mare Deisenroth. Data-efficient reinforcement learning with probabilistic model predictive control. In Proceedings of the Twenty-Fiin First International Conference Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1701-1710. PMLR, 09-11 Apr 2018
[32] James Kelliher. Stream functions for divergence-free vector fields. Quarterly of Applied Mathematics, 1914)1163-174, 2021.
[33] Jonas Kiessling, Emanuel Ström, and Raúl Tempone. Wind field reconstruction with adaptive random fourier features. Proceedings of the Royal Society A, 477(2255):20210236 2021.
[34] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6988. 2014.
[35] Markus Langge-Hegerman
Algorithmic linearly constrained Gaussian processes. Advances in Neural Information Processing Systems, 31, 2018.
[36] Markus Lange-Hegerman Linearly constrained Gaussian processes with boundary conditions.
In International Conference on Artificial Intelligence and Statistics, pages 1090-1098, PMLR. 2021.
[37] Markus Hegerman and Daniel Robertz. On boundary conditions parametrized by analytic functions. In International Workshop on Computer Algebra Scientific Computing, pages 225-245. Springer, 2022.
12

---

[38] Yann LeCun, Bernhard Boser, John Denker, Donnie Henderson, Richard Howard, Wayne Hubbard, and Lawrence Jackel. Handwritten digit recognition with back-propagation network. Advances in neural information processing systems, 2, 1989.
[39] Mario Livio. Why symmetry matters. Nature, 490(7421):442 473, 2012.
[40] Lennart Ljung. On consistency and identifiability. In Stochastic Systems: Modeling, Identifica- tion and Optimization, I, pages 169-190. Springer, 2009.
[41] Svenja Lowitzsch. Approximation and interpolation employing divergence-free radial basis functions with applications. Texas A&M University, 2002.
[42] Xiaoyu Lu, Alexis Boukouvalas, and James Hensman. Additive Gaussian processes revisited.
In International Conference on Machine Learning, pages 14358-14388 PMLR, 2022.
[43] Jun Ma, Fuqiang Wu, Wuyin Jin, Ping Zhou, and Tasawar Hayat. Calculation of hamilton energy and control of dynamical systems with different types of attractors. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(5), 2017.
[44] Andrew J. Majda and Andrea L. Bertozzi. Vorticity and Incompressible Flow. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2001.
[45] Filippo Maria Denaro. On the application of the Helmholtz-Hodge decomposition in projection methods for incompressible flows with general boundary conditions. International Journal for Numerical Methods in Fluids, 43(1))43369, 2003.
[46] Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman. Invariant and equivariant graph networks. In International Conference on Learning Representations, 2019.
[47] Kaspar Martens. Enabling feature-level interpretability in non-linear latent variable models: a synthesis of statistical and machine learning techniques. PhD thesis, University of Oxford, 2019.
[48] Takashi Matsubara and Takaharu Yaguchi. FINDE: Neural differential equations for finding and preserving invariant quantities. In The Eleventh International Conference on Learning Representations, 2023.
[49] Takashi Matsumoto. A chaotic attractor from Chua's circuit. IEEE Transactions on Circuits and Systems, 31(12):1055-100 1058, 1984.
[50] Francis J Narcowich and Joseph D Ward. Generalized hermite interpolation via matrix-valued conditionally positive definite functions. Mathematics of Computation, (3(208):6116687 1994.
[51] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, tegforeeformance deep learning library. Advances in neural information processing systems, 32, 2019.
[52] Jim O Ramsay, Giles Hooker, David Campbell, and Jiguo Cao. Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 6955:74117796, 2007.
[53] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 11 2005.
[54] Katharina Rath, Christopher Albert, Bernd Bischl, and Udo von Toussaint. Symplectic Gaussian process regression of maps in Hamiltonian systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(5):03312,, 2021.
[55] Marco Reisert. Group integration techniques in pattern analysis: a kernel view. PhD thesis, Freiburg (Breisgau), Univ., Diss., 2008
[56] Marco Reisert and Hans Burkhardt. Learning equivariant functions with matrix valued kernels. Journal of Machine Learning Research, 8(15):385-40) 00885-408,40 2007.
13

---

[57] Jack Kichter-Poweld Yaron Lipman, and Ricky T. Q. Chen. Neural conservation laws: A divergence-free perspective. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022.
[58] Steffen Ridderbusch, Christian Offen, Sina Obbr-Blöbaum, and Paul Goulart. Learning ODE models with qualitative structure using Gaussian processes. In 2021 60th IEEE Conference on Decision and Control (CDC), pages 2896-2888 2896. IEEE, 2021.
[59] Magnus Ross and Markus Heinonen. Learning energy conserving dynamics efficiently with Hamiltonian Gaussian Processes. Transactions on Machine Learning Research, 2023.
[60] C Sarasola, FJ Torrealdea, A d' Anjou, A Moujahid, and M Graña. Energy balance in feedback synchronization of chaotic systems. Physical Review E, 69(1):01 1606, 2004.
[61] Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E(n) equivariant graph neural networks. In International conference machine learning, pages
PMLR, 2021.
[62] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inpure In Y. Weiss, B. Schölkopf, and J. Platt, editors, Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
[63] Ercan Solak, Roderick Murray-Smine WE Leithead, D Leith, and Carl Rasmussen. Derivative observations in Gaussian process models of dynamic systems. Advances in neural information processing systems, 15, 2002.
[64] Steven H Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press, 2018.
[65] Tomoharu Suda. Construction of Iyapunov functioos using helmboltz-hooge decomposition. Discrete and Contimuous Dynamical Systems, 39(5):2477-2454, 2019.
[66] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
[67] Peter Toth, Danilo Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins. Hamiltonian generative networks. In International Conference on Learning Representations, 2019.
[68] lari Vallivaara, Janne Haverinen, Anssi Kemppainen, and Juha Röning. Simultaneous local- ization and mapping using ambient magnetic field. In 2010 IEEE Conference on Multisensor Fusion and Integration, pages 14-19, 2010.
[69] Ilai Vallivaara, Janne Haverinen, Anssi Kemppainen, and Juha Röning. Magnetic feld-based SLAM method for solving the localization problem in mobile robot floor-cleaning task. In 201 15th International Conference on Advanced Robotics (ICAR), pages 198-203, 2011.
[70] Pantelis R Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos. Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics. Neural Networks, 22:1191-117 2020.
[71] Niklas Wahlström. Modeling of magnetic fields and extended objects for localization applica- tions. PhD thesis, Linköping University Electronic Press, 2015.
[72] Niklas Wahlström, Manon Kok, Thomas B. Schön, and Fredrik Gustafsson. Modeling magnetic fields using Gaussian processes. In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 3522-32220003 222-3526,222 2013.
[73] Nils Wandel, Michael Weinmann, and Reinhard Klein. Learning incompressible fluid dynam- ics from scratch- towards fast, differenntable fluid models that generalize. In International Conference on Learning Representations, 2021.
(74] Maurice Weiler and Gabriele Cesa. General E(2)-equivariant steerable CNNs. Advances in neural information processing systems, 32, 2019.
14

---

[75] Holger Wendland. Divergence-free free kernel methods for approximating the stokes problem. SIAM Journal on Numerical Analysis, 47(4):3158-3179 2009.
[76] Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapo- lation. In International conference machine learning, pages 1067-1075. PMLR, 2013.
[77] Ping Zhou, Xikui Hu, Zhigang Zhu, and Jun Ma. What is the most suitable Lyapunov function? Chaos, Solitons Fractals, 150:1 11144, 2021.
[78] Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence. Kernels for vector-valued functions: review, Foundations and Trends@ in Machine Learning, 4(3):195-266, 2012.
15

---

Appendix
A HHD and dynamical systems
This section provides brief introduction to the connections between the Helmholtz-Hodge Hodge decompo- sition (HHD) and some dynamical systems needed to understand this work.
A.1 HHD and dissipative Hamiltonian systems
We begin with a brief review of Hamiltonian mechanics For dynamical system with degrees of freedom, the Hamiltonian formalism describes the system by defining scalar function H (x) known the Hamiltonian, where the system state x (q. E IR2N is described by generalized coordinates E IN and E I N in the phase space, corresponding to generalized position and
momentum, respectively. The time evolution (q. p) governed by the symplectic gradient of its 8H If the Hamiltonian is not explicitly time-dependent, [4, p!
Hamiltonian, i.e., q 3p
0H P
3q
ah defines a divergence-free vector field oq Öp 0p å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, account for energy dissipation, term responsible for the volume contraction
the phase space can be introduced into the Hamiltonian dynamics. This term represents the energy lost by the system due various dissipative forces and is typically modeled by gradient fields of some scalar functions (-), so the motion equations of the dissipative Hamiltonian system be given by
TH
OH
f (q.
(q())
(19)
0p
3q_
curl-free
divergence-
In dissipative Hamiltonian system, the Hamiltonian still represents the system's 's total energy, but the additional damping term causes the system to lose energy over time. One common example of Qp. where Q E NNN is
the damping term is induced by the Rayleigh function, fincctin,,, D (p) symmetric positive-definite matrix called the Rayleigh dissipation matrix. The Rayleigh function friction, air resistance, and viscos- ity-in the context of Hamiltonian mechanics. When the Hamiltonian H (.) represents mass instead of energy, Eq. [9 can also be used to model 2-dimensional compressible fluid dynamics, such the dynamics of ocean current fields.
provides elegant way include dissipative forces-sude
A.2 HHD and chaotic systems
The dissipative Hamiltonian system (Eq. 19 is in an explicit form HHD, because the Hamiltonian vector field is divergence-free and the dissipative field V D is curl-free. However, as its divergence- free part is governed by the Hamiltonian equations, Eq. 19 can only used to describe a subset of even-dimennional systems. To extend its scope to study chaotics systems, (60) proposed to approach energy for dimensionless dynamical systems by using HHD as follows,
x = (x) = fase (x) + feart (x), x E R"; V H..... (x) = 0; H TH (x)
(20)
where faie and feurt are divergence-frn free and curl-free components of a dynamical system, respectively. According this criterion (Eq, 20)). HHDs of chaotic systems were approached calculate the energy for synchronizing two chaotic systems [60], to analyze the stability of chaotic systems (77).
and to design energy modulation based controllers 4
B GPs and linear operators
Gaussian processes (GPs) are closed under linear transformations (ef. Lemma 2.2 [36] Lemma 2.1 in (25).. Let L be linear operator acting realizations of GP (Mg (x), Kig (x,x)), then under
16

---

the assumption that 2 commutes with expectation, the mean of Lg is given by
IE [Cg (x) LR. (x)) = Lig (x),
(21)
and the covariance is such that cov [Gxg (x), fx'g (x')= (Ca9 (x) Exl1 (x)) (fxe (x') - fxe (x))
Cxl (9g (x(x) - 1 (x)) (g
- - -x(xx)) cx
(22)
= Cxtig (x, x') Cgp.
A common application of this technique is construct GPs with realizations in the solution set of linear differential equations, assuming that C is a linear differential operator [29] 36, 7 25]. And similarly, 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. has been widely used, its validity as covariance function is rarely discussed in existing works. In GP framework, a kernel x (x, (x)x') RR x IR + R is valid covariance function if is:
(i) Symmetrie, i.e., (x, x') (x', x) Vx, x' IR", and
(ii) Positive semidefinite, i.e., Eij cf (x;, Xj) Cy 2 O, for any finite set {x:} C I and {c:) C R
These two conditions for kernel can be verified if and only if there exists feature map Φ (x) such that & (x, x') (x) Φ (x') 53 So, if Kig (x, x') is valid kernel, it holds that
Cx (x. CT = Lxx (x) (x') 4 (CxØ (x)) (x)))
(23)
Therefore, the transformed kernel Lx (x,x') LT is guaranteed to be valid covariance function provided that its underlying kernel fig (x, x') is. However, when L is a differential operator, Kg should be twice differentiable, which is satisfied by most of the standard kernels, such as the squared exponential kernel.
C Basics for Euclidean group
This section gives the basic definitions about the Euclidean group. In the context of this work, the most important example a symmetry group the Euclidean group E (n) or its subgroups. The set of all elements in E (n) can be denoted
E (n) = {(A,b) (n) ER
(24)
where (n) = {AE R | AA = } is the orthogonal group. Any element g = (A, b) E (n) represents a translation followed by orthogonal transformation, the action of g point R" is given by a linear mapping:
Ly
R x x
b).
(25)
Therefore, E (n) comprises all isometries of a Euclidean space, i.e. for all x, x' E IR" and E E (n). have
Ly (x) - Lg (x')l -x - x1||
(26)
where I1-I1 is the Euclidean norm. All intuitive geometric transformations in R" can be described by subgroups of E (It), such as
1. Translation: The group of all translations in R" is denoted by (R" +). For any E R",A translation a transformation that moves a point R" by Ly (x) =v
2. Rotation: The group of all rotations in IR" is represented by the set of special orthogonal matrices SO (n) = R E O (n) det. R 1}. where rotation matrix R transforms point by T.n (x) Rx.
3. Reflection: Reflections in Rn forms subgroups of the orthogonal group O (n). Reflections correspond to mirror symmetries. They mirror points across hyperplane. For a hyperplane with a unit normal vector n, the action of reflection is defined as Ln (x) 2(x n)n.
17

---

D Proofs
This section restates and proves the theorems in Section 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 be Euclidean group or its subgroup, and let V IR - R be G-invariant scalar function. Then, the curl-free vector field feart : IR + R2" defined by feart (x) = VV (x) is G-equivariant.
Proof. To prove that VV (x) is G-equivariant, by Definition 5.1 we need to show that for any x E IR" and y G. we have
VV Lg = J2, VV,
(27)
where Lg (x) : R > R" is the action on x, and Jr, is the Jacobian matrix Lg. To show this, first recall the chain rule of the gradient operator,
(V Lg) IL (VV Lg),
(28)
where JT denotes the transpose of L,, Since V (x) is G-invariant, for all 9 we have
== Lg.
(29)
Now taking the gradient of both sides, and combining with the chain rule, we obtain
VV =V= V(VoLV- J2, (VV o Lg).
(30)
Then by multiplying Jr, on both sides, we have
JA VV J., (VV Lg)
(31)
Considering G is a subgroup of the Euclidean group, the action of any group element E 9 on IR can be represented by
Lg (x) = (x
(32)
where A E O (n) is an orthogonal matrix. Therefore, Jr., J = AAT = 1, where / is an identity matrix, so we have J,,, VV = VV Lg. which means that the curl-free vector field fourt = VV is Geequivariant, as desired.
D.2 Proof of Theorem
5.3
Theorem 5.3, Let G be a Euclidean group or its subgroup, and let v ; IR" - Rn be a G-equivariant vector field. Then the divergence-free vector field faio = (V A....., V A, is G-equivariant, where A; denotes the i-th row of the skew-symmetric matrix-valued function A Jv J
Proof. To prove this theorem, we need to show that the divergence-free vector field faiv satisfies the geequivariance condition, ie., for all E G and x E R
fiiv
Lg JL, faive
(33)
where Lg (x) R² + Ra is the action of g on x, and r.,, is the Jacobian matrix of Lg-
Define the divergence- free vector field fair V A1,..., V AT where A; is the ith of the skew-symmetric matrix field A = J, JJ and Jv is the Jacobian of some vector field v R R Given that v G-equivariant, for all g E 9 and x E R we have
voLg V.
(34)
Now computing the Jacobian of both sides, by the chain rule of the Jacobian, we obtain
JvoL, = (Jv o Lg) = Jv.
(35)
Then applying J[ both sides, we have JL, (Jv Lg) Jr, J J Jv. Since G is a subgroup of the Euclidean group, the Jacobian of the action of any group element g G on R2" is orthogonal
18

---

matrix, i.e., J1, E EO (n). Therefore, JL JLg = I. where I is an identity matrix, so it holds that
Jv -1, (Jv bg) lig Then by substituting into Jv - J we obtain
Jv -J. = II, (Jv l J L.g) Jr., Jh, (A Igg) where the equation at the (i,j)-th entries of both sides is given by
(36)
a (Lg)e Σ a (Lg))
Aij
(AM L.g) Ōrg
(37)
Or:
k=1
1=1
Then by substituting Aij into the construction of fain. for any E G and x E IR", we have
fdi
V
A1
A
8A1 dz,
Σ Ej=1
aAni Ozy
Ej=1 Ek=1
8(Lg)
0z1
-i=1 day
o(Lg) (g)
d2y
Σ j=1
8(Tg) ((AuoLo)
Oin
1==1
Biy
day
8(1gg
((1gh
((A
2=z=1 Ei=l
k=
dey
day
day
Σ 1
0(L.)
2f==1 Σ
a(Lgh ((ggoo 0zj
k=
02n
Omg
Σ
((4g)
2j==1 E
"
0(Lg)
Σ n p=1
0(La),
(An o Lgg dep
0z
ow;
(38)
Σ
8(Lg)e
2j=1 E1=1
8(Lg)
Σp=1
0(Lg)
AAn o Lug
Ctp
Σ k=
K(Lg)
n Ep=1
aAx o 1ug) S2,
Σド=1
d(bg) a(Lo) 02, 02,
021
Ek=1
8(Lg)e
Own
21=1 Ep=1
Tig)
Σラ=
8(Lg) o(Lo)
Ozs 0eg
E
8(Lu))
Ei=1
(AAli o Log)
Ek=l
a(La)k
Ei=1
AAkI
Lgg Lg)]
di
Ozn
Ozl
E
8(La)k
(V- Aa) o Lgg
Ek
a(aa) (V AA)
dm
dan
JT, [(V A1) o by
AV- Lg]
=J[, (faio I.g).
Then, applying Jr., on both sides, we obtain Jr., faiv = faiio T.g, which completes the proof.
E
Physical systems in the experiments
This section reviews the physical systems used in our experiments, listing their governing equations, HHDDs and symmetry properties.
E.1
Dissipative Hamiltonian systems
The connection between HHD and dissipative Hamiltonian systems is reviewed in Appendix A. In experiments, the models evaluated with the following two dissipative Hamiltonian systems.
Damped Mass Spring A damped mass- mass--pring system is mass attached to spring that oscillates periodically around an equilibrium position. Its Hamiltonian in natural units is given by
H (g.p)
(q² + +2²)
(39)
19

---

where q IR is its position, and is the momentum conjugate to The Hamiltonian represents the total energy of the oscillator. Without energy dissipation, the motion of the oscillator described by the Hamiltonian vector field (div-free vector field face) derived from the Hamiltonian:
Edii
1 P
HH
p
(40)
0q
Then we show that the Hamiltonian vector field has SO())-qquivariance Apply a rotation by an angle 0 that transforms any state vector [g. p] to p''' where o' = cos0-ppsi 0, p' qgsin 0 cos 0, then for all 0 E R, the equivariance condition (Eq. 14] can be easily verified by
1 p' - L-
cos 0 sin 0 sin
coos 0 sin 0 sin 0 cos favo (g.p).
Easo (q' .3)
p
(41)
cos
Therefore, the Hamiltonian vector field is equivariant with respect 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 fowrt) is induced
by the Rayleigh function D (p) = =pp2/2 /2:
fourt = V (g.p) D = [0,-pp]"
(42)
where p is the friction coefficient and set it at 0.1. This dissipative field has the symmetries of the g-axis translation (feuri (q + 9: p) feurt (q. p). 4g E R) and the rotation by an angle of π radians (the so-called odd symmetry, i.c., fewrt 9- vector field (Eq, 42) and the Hamiltonian vector field (Eq. 40). the dynamics of the damped harmonic
f (g. p)). Then by summing the curl-free
oscillator is characterized by
= faiv + feurt - [p, -q 0.5p)
(43)
which only demonstrates the odd symmetry.
Damped Pendulum A damped pendulum is physical system consisting of weight suspended from a pivot, subjected a resistive force that gradually reduces its oscillation over time. By defining its state in the phase space through [g. p] E IR? the dynamics of a damped pendulum and its HHD can be given by
8H
an
3p
DH
3g 2mgl sin q
P
+
=
m mmglan sin pp]
(44)
PP
0q
8p
divergence- free
curl-free
where the Hamiltonian H and the Rayleigh function D are given by
12p²
1
(g. p) = 2mgl (1 cos q) + 2m (p))
(45)
In experiments, the gravitational constant g, the mass m, and the pendulum length / were set as g = and m = / = 1, and the friction coefficient was set as p = 0.1. The dynamics of the damped pendulum, well 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.
E.2 Chaotic system
The Chua circuit 49 is simple electronic circuit that exhibits chaotic behavior, and it has applica- tions in various fields, such as secure communication systems and random number generators. The ODE of a Chua circuit with its HHD is given by
α (y x3 cz))
azy T +
a
( (-23
C.
x - y +
(46)
By
3y
divergence free
cuel-free free
where x, y, 2 are the phase space variables, and a, e, 3 are the system parameters, which were set to 10, 14, and 16 respectively in the experiments. The Chua circuit is equivariant under
20

---

rotation about the origin: Lx (.. + (-r.-y,-2) If the Chua circuit is decomposed by the HHD in Eq. _o more symmetries are exhibited. Specifically, its div-free component inherits the odd symmetry, but additionally presents the translation invariance along x = The curl-free vector field
invariant along the z-axis and equivariant under mirror reflections across the coordinate planes,
L (I, y, 2) + ((-1)' x, (-1)' y, Vi,j E (0,1)-
(47)
This HHD not only uncovers more knowledge of the Chua circuit's 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:
0H
0H dy
0H
ay dr
+ (x+ + 2)
By dz = 0,
(48)
which is satisfied by the quadratic form:
1
1 B
H
- x y
(49)
2
a
As presented by (77), this energy function indicates that the Chua circuit keeps oscillatory when the energy release along the æ-axis is enough to balance the energy pumping along the y-axis and z-axis.
F Evaluation metrics for learning ODEs
RMSE and MNLL of state derivatives To evaluate the models' performance in terms of learning = f (x), we first compute the woot mean squared error (RMSE) between the predicted time derivatives &; and the ground truth x,, i.c.
xIIe) $
m
RMSE
x:
(50)
171
where ||-112 is the Euclidean norm, and is the number of test data. The lower the RMSE, the better. The test set [(x,, xt, } was generated by sampling grid of points with the resolution of 10. In addition to the RMSE for evaluating the regression results, further use the mean negative log likelihood (MNLL) to evaluate the prediction uncertainty provided by the GP models, and is calculated by
m
MNLL
Σ
log, Nox: XR, var (x:)),
(51)
m.
i=1
where var (x:) is the prediction variance at a test state x,. 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 770 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 at which the normalized root mean square error (NRMSE) between the predicted state &; and the ground truth exceeds given threshold e,
1
VPT T arg min (NRMSE (&,, x1) GYDSEST),
(52)
where NRMSE x:) (x) -(x, 1 E diag (1/01,..., 1/o,), with o
x:)
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 trajectory from an initial condition was solved by the Dormand-Prince method (dopri5) [[S] implemented in torchdiffed" integrating forward in time at a frequency of 25 Hz. for 15 seconds, with the relative and absolute tolerances of 6
21

---

G Implementation details
The experiments were performed on single Nvidia GeForce GTX 3050 Ti GPU, and all of the models were implemented with PyTorch [51]. The GP-based models (Ind-GH GIM-GP, Div-GP, HHD-GP and SPHHD-GP) were trained by maximizing the log of their marginal likelihood:
- Y (K + E)+
logp (Y X) log N (Y
log K + Σl log 2n,
2
(53)
where |-| computes the determinant of the covariance matrix + Σ. Training 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 with learning rate 0.01 for 3000 gradient steps.
34
G.1
Ind-GP
The Ind-GP models each dimension of a n--diennsional dynamical system independently, so its matrix-valued kernel is constructed by
Kind = diag Ka))
(54)
where Ni, i 1,..., independent scalar kernels. And standard choice for Hi is the squared exponential (SE) kernel,
Uige (x,x)= = 02cxp (-11-211
|| x(1²)
(55)
which has two parameters: determines the variation of function values from their mean, and / controls the length scale on which the function varies. Realizations GPs with SE kernels are dense in the set of smooth functions O (R" IR) (cf. Prop. in (37)) For fair comparison, also used the SE kernel (Eq, (55) to construct kernels for the other GP-based models.
G.2 GIM-GP
The GIM-GP produces predictions with the desired symmetry. Each of the systems in experiments exhibits the odd symmetry as a whole. So their symmetry groups can be given by Godd = f - where is the n-dimensional identity matrix, and the group elements are linear representations of the group actions, i.e., Lg (x) gx, Vg E G. Therefore, according to Eq. 17 the GIM-kernel for Godd constructed by
Kigim (Nae (x,x') (x,-x') -x')) Inx
(56)
where the SE kernel (Eq. 55 was used as the basis kernel x: in Eq. 17
G.3 Div-GP
GPs with a div-free kernel can be used to approximate conservative dynamics. According to Eq, two-dimensional div-free kernel was constructed by
12
8 8 (x,x²)
x (xx,x' Cungchar
Edio
82 (x,x')
89 Ne (x,xx
(57)
Dury Darl
Chan Ourl
where Ku in Eq. 12 is denoted by Kи 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. used for fi.. The partial derivatives 57. [57werr were calculated by automatic differentiation in PoToch,, so did not need to derive its analytic expression.
Then we constructed the div-free kernel for the Chua circuit system. According to Eq.
three-
dimensional div-free vector field is given by fav (x) Vu (x), where v [112 113 23 E 33x3
22

---

with its components given by
0
0
1/12
Or
Ōz
(58)
0
0
0
Ong
0
0
17
ORs
13
Oay
(59)
-
0
0
0
0
0
2/23
0
0
Ozs
(60)
0
-1
0
Then by assuming = [u12, 413, U23 for the Chua circuit according to Eq.
GP (0, Ku
Kind), we constructed the div-free kernel
12
Ozz
Ozs
61
0
O2z
.z1
Ridio
0
8 23
0
62
0
0
(61)
Orr
0
N3
0
Oz
80 (xxx²)
)
+ Pe (xxx')
@rra Sera((x (x.x')
Sr (x.xx'
8 11 (xxxx')
02 (e((( (x.x²)
@
(xxx'
o Pnaa (xxxx')
(62)
des der
82 82 (x.x))
8 0P (x.x')
8 Pma (xxxx'))
8f (xxxx)
+ Cerathigg
dzz
duy
where K1, and Ky are all independent SE kernels (Eq. 55
G.4 HHD-GP
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 build curl-free kernels according to Eq. X For the damped mass spring and the damped pendulum, their two-dimensional curl-free kernel was constructed by
82 xxxv((xxx)
8uy (x,x') Cxrl puv (x.x') Ooiguall
Ozitha 8" ('nv(xx ('vv(x.') 0x2dhs
Beurr
(63)
Similarly, the three-dimensional curl-free kernel for the Chua circuit was constructed by
@ev (x,x') Ou, Oue d' (x.xx)
@ Øerv (xxx))x)) Ox 0x,
@2 (x,x))
Oxx dury @xv (xxx)
0 kyy (x,x²)
Keurl
(64)
8ev (xxxx)
@ By (x,x)
Oxy (xx,2))
where By is SE kerne..
G.5 SPHHD-GP
G.5.1
Symmetry--peserving curl-free kernels
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, G-equivariant curl-free kernel is obtained by constructing its potential kernel Ky according to [6
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, (x) x + (g. 0), for all E R; (ii) rotation by angle of radians, i.c., Lg (x) gx, for all g E (12, - -12)
23

---

These symmetry groups were enforced to a SE kemee sequentially according to Eq. 16 where [g can be ignored. Specifically, the translation invariance was enforced using the Gaussian integral formula:
KV =
Kse
(x (..'+(0.0') dg
(65)
/ exp( -- 9)²
Kse (p.p')
)gg
(66)
212
= =V2llssee (p.P),
(67)
then based on which the и rotation invariance was enforced by
Kv /2Tl (Nee (p,p') Kise (p, -p')).
(68)
The curl-free vector field of the Chua circuit (Eq. 46) has similar two types of symmetry: (i) translation along the z-axis, i.e., Lg (x) x + (0, 0, g), for all E IR;
(ii) mirror reflections across the coordinate planes, i.e., Lis (x, + (1-1) x, (-1)² y, 2) for all i.j {0, 1).
Therefore, By for the Chua circuit was constructed by Keze (x, gx') G diag ((-1), & e (0.1)).
Kv V /Tm
(69)
SEG
G.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 with the same equivariance. For the div-free vector field of the damped mass-spring system with SO(()-equivariance, its Ky was constructed by the GIM-kernel (Eq. 17
coss sin A
sin A cos 0
Ki (x,
Kise (x, g6x') god0, where ge
(70)
Then this Ky (Eq.
70 is substituted into Eq.
18 to construct KH for its div-free kernel Kaiv (Eq.
57
0²
0²
0² [kv|1,2
02
KH Oxzdnz
[Kv]1,1
[kv12,2
0x, [nvl2,1
(71)
0r,0x
Oxzdzl
212 ||x gox'll
10 (x, g6x')
do
(72)
which admits closed-form solution, by sampling discrete rotations of { in n = O, o....,7! 7]. group. And by setting / (0, m Eq. 72 was used to construct кı for the div-free vector field the damped pendulum, which has the equivariance under rotation (odd symmetry).
used numerical approximation of the integral (Eq. 72
whose rotation matrices form a finite
The div-free vector field of the Chui circuit (Eq. 46) has the odd symmetry and the translation symmetry along = z. Ky respecting the translation symmetry was constructed by Ky = K 13. where K is given by
K fige (e, L1,, + 1g, -91-) dg
0..0,-991
(73)
Koc (y,1²) 1::0 exp(
(x - x' 9)2 + (z z'' + 9)² dg
(74)
212
to0
2 g
mt - x - Bg
(m)
2 t2g
2
= Kac (y,0²)
exp
dg
(75)
212
/ too
(- (2) - o 23 xg/)²
= =Vilisse 2b)
exp
dg
(76)
412
Villnsse (p.p').
(77)
24

---

where p [/2 (x ( 2) And further enforced the odd symmetry by
R val (Rise (P,P') Kisc (P.-P')
(78)
Then, Ky Kz was substituted into Eq. 18 to construct Nin by
=0x(xxx)
(xx(xx')
8x(xxx')
o ( Bard o² xxx(xxx) Harittrr & (x((x.)))
Oagiday
dritzi 0 0m((xx'' d2g 0 xxx(xxx')
Dazday
0 (xxxx')
xm(xxx))
Ku
(79)
Ongttag (x((x.'))
021dz
o xx((xxx)) OzzdEz
+
Finally, this Ku (Eq. 79 was substituted into Eq. 12 to construct the symmerry-preserving div-free kernel for the Chua circuit.
G.6 D-HNN
We used the released code 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 KH for constructing their div-free kernels (Fq.. (5)) can be interpreted as placing GP prior the Hamiltonian function. Therefore, joint GP describing both the Hamiltonian function (x) and the dynamical system (x) is given by
H (x) [(x)]
(!)
KH (x, x') RiGH (x,x)
KH,F x Kihd (x, x
GP
(80)
where KH. (x,x) Ne, H x) with
KH,F (x, x') = cov [H (x) (x)) ] = 10p -0,-] KM (x, x')x;);
(81)
81 (x,x') x') cov f (x) H (x') 10p, -0,] KH (x,x').
(82)
After training the model f (xx GP (O, Kand (x,x') x')) on noisy observations Y Y (y1, (y.....,, ym at
states X = [xi,..., Xm we are interested in predicting the value of Hamiltonian function H (x.)
test state x.. Since these data determine (.) only up an additive constant,
that we have an anchor point H (Xo) which can be chosen arbitrarily. Then, according to the GP prior
(Eq, 80). H (x.) and Y и [H (No) YJ jointly distributed as
(x.) Y
~ N (1o).
KH (X,, x.)
k K
(83)
kT
where k = [N (X.; Xo): KHH (x., X)], and
KH (Xo, Xo) Ki. (X, xo)
KH.P (Xo, X)
K
(84)
Kand (X, + tl
Then, we obtain the posterior distribution
p (H (x.) Yx) = - (kK- YYn KH (xe, (x),Xx)= ) - kK kk-
(85)
where the mean function is used for energy prediction.
y, f f(x)) c,c 25 (0,
25

---

G.7.2 Predicting energy of the Chua circuit
60] stated that the div-free component of a chaotic system is always orthogonal to the gradient of its energy function, i.c., V fdiv (x) = 0, for all x E R The energy function of the Chua circuit in quadratic form:
H = (11-1100 ---- 2 a
(86)
B
So we parameterize the energy function of the Chua circuit by H (x) = 2 + a2y² + az2
where (a), 022 03] are parameters. Then by learning a divergence-free vector field faso (+) 2
through our model, we can estimate the parameters by minimizing Σ. (V)TT( (
at finite
number of sample points { (m = 500 in our experiments). And to climinate the solution at [0, 0, 0j. we add an equality constraint at random point of the ground truth {xo, H (xo)). Therefore, the parameters of H are solved in a convex quadratic program (QP):
min a Qa
(87)
s.t. H (xx) H (xu)
where Q - q qoq with q [aai (x.),....ff.. Eain (xm)
E mRmx3, And in the experiments, we solved
the QP using the solver provided by CVXOPT
H Uniqueness and Symmetries
The two constituent kernels of our model define the space of divergence-free vector fields (faiv E Fase) and the space of curl-free vector fields (feurt E Feurt), respectively. These two spaces overlap partially due to the presence of harmonic vector fields (fharm E Faiv n Fourt) 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 Gdiv and Geurt. Therefore, the uniqueness property of HHD depends on the space of harmonic vector fields that respects the union of two symmetry groups. i.c.
(88) If F Fharm= 0, the HHD is unique. i.e. our model is identifiable. A harmonic vector field can always be represented by the gradient field of harmonic function h, which is a scalar function satisfying V Vh- = 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 rotation group Gáse SO (2) and translation group Gourt = {(g, 0) g E RJ. A harmonic vector field frorm. V that satisfies this translation symmetry implies that its harmonic function h (q.p) is independent of the variable g. so the harmonic vector field is given by frarm = (0, 0,h) and 0,h should be constant to satisfy the Laplace's 's equation, but the harmonic vector field in the form of constant clearly contradicts rotation symmetry Gávv. Therefore, there is no harmonic vector field that respects both the symmetry groups Gaiv and Geurt. Similar conclusions damped pendulum and the Chua circuit. The Laplace's '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.
Fherm = {fharm fharm o In = fharm, VX RR" Gdiv U Geart)
be drawn for the
I Computational complexity
Compared to the diagonal kernel that independently models each dimension of a dynamical sys- tem 1 3T 68 69, our kernel can capture correlations among output dimensions because its off- off-diagonal elements are non-zero. However, this advantage comes together with a low computa- tional efficiency. The exact inference of HHD-GP over an n-dimensional system with m data points has cost of O (m²n²) in time, due to the inversion of its full covariance matrix K, E Rmmxmm 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 [22 66), reducing the computational complexity to O ( [M m where M (M < m) the number of inducing point to approximate the covariance matrix.
26

---

J Additional Experimental Results
J.1 Results with increasing number of training data
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. []] As expected, the performance the models improves the amount of training data increases, but can observe that our model (SPHHD-GPP performs best at every data amount relative to the baseline models, verifying the robustness of our model against data sparsity.
Mass-spring System
Mass-spring System
Mass-spring System
100
HDG
10
RMSE
VP
MNLL
10-2
10-1
Number data
Number
Pendulum
Pendulum
Pendulum
Ducc
10a
ki
Ho-2F bb:2
10
r&
RMSE
VPT
MNLL
10-1
100
200
300
400
100
200
300
400
100
200
300
400
Number of data
Number of data
Number of data
Figure 5: Results with an increasing number of training data (20, 100, 260, 420), which are averaged over 10 independent experiments performed by resampling the training sets and model initial param- eters. The first and second columns are the results for the mass-spring system and the pendulum, respectively.
J.2 Trajectory prediction of the mass-spring system and the pendulum
Fig. Gand Fig. present the trajectory predictions for the damped mass-spring system and the damped pendulum, respectively.
J.3 Example Decomposition
Fig. 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. M(a)). they learn completely different decompositions (the second and third columns in Fig. ((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. [bb] 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.
J.4 Experiments of predicting dynamics with unseen frictions
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. the curl-free dynamies 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 p = 0.1 (p is the friction coefficients), then generalizing the trained models to dynamics when p = 0.05 and p = 0.5 by multiplying the learned curl-free component with constants 0.5 and 5.0, respectively. The results
27

---

1.0
0.15
0.10
IHD
0.5
4
die4 op
199 oP
-0.05
-0.5
-0.10
20
30
-0.15
40
time
1.0
0.15
0.10
0.5
p
Oround Tuin
Oround Truth
0.5
Di--p
-0.05
-0.10
-0.15
time
Figure 6: Comparison predicted trajectories of the mass-spring system.
0.75
0.50
0.25
9
Binund wurh Hnm o7 Saiin
0.00
Grand Tran
IHE-67
0.25
Dirkeu Indice
0.50
-0.75
10
15
20
19
20
21
22
time
time
2-
isa Gai--b--i di Ca--ca
Hep--p
spe--t-s Ba-ch Denk
-2
CIM-C2
CiNicr
20
100
20
21
22
time
Figure Comparison of predicted trajectories of the pendulum.
shown in Fig. 9 where observe that the three models accurately predict the trajectories of the systems with a friction coefficient 0.1 (first column in Fig. 9 but HHD-GP and D-HNN are difficult generalize cases with friction coefficients of 0.05 and 0.5 (second and third columns in Fig. 99. In contrast, SPHHD-GP effectively captures the dynamics of different friction conditions, even though 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.
J.5 Results of increasing noise level
To evaluate the robustness of the models, added the experiments of increasing the standard deviation {0.01, 0.05, 0. 0.10, 0.20) of Gaussian noise in the training data. Fig. proseetstss the results on the pendulum. As expected, the performance of the models degrades as noise increases. However, can be observed that SPHHD-GI still performs best at all noise levels compared to the baselines, meaning that our model more robust data noise. Furthermore, from Fig. [O(d] we observe that the energy predictions of HHD-GP and D-HNN generate consistently significant errors, insensitive to noise level 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
28

---

HHD-GP GP
3.5
2.7235
2.7230
3.0
5.7985
2.5
2.7275
2.0
SPHHD-GR GP
1.5
1.0
(a)
(b)
Figure 8: Example predicted HHD of the damped mass-spring system by HHDD-G (the first row) and SPHHD-GP (the second row). (a): predicted vector field (the first column) with its div-free (the second column) and curl-free (the third column) components; (b): the associated variance of the div-free (the first column) and curl-free (the second column) predictions.
p= 0.1
0.05
p= 0.5
Greadd
Mass Spring
des Fur
Full
21100
HIOLD
Pendulum
Grund Pus Iac- Da
Crond Push IIO
p
AM
220
41
420
Figure 9: Adapting the trained models to predict the trajectories for different friction coefficients. The models are trained under / = 0.1, where p is the friction coefficients, then they are generalized to predict trajectories for م 0.05 and with constants 0.5 and 5.0, respectively. The lines the figure represent predicted trajectories from different models (distinguished by different line styles), with the initial points for the two systems (the first and second rows) being (1, 1) and (1.5, 1.5), respectively. The background vector fields are the ground truth dynamics.
0.5 by multiplying the learned curl-frec component
noise, further demonstrating the effectiveness of addressing model identifiability through symmetry constraints.
J.6 Trajectory and energy prediction of the Chua circuit
Fig. 12 presents the comparison of trajectory predictions for the Chua cireuit, and Fig. 11 compares the energy prediction along the trajectory.
29

---

RMSE
MNLL
f
H-NN
10
F
GIM-GP
f
HHD GP
-2
Indi
F
SPHHD-GP
-3
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
SD of noise
SD of noise
(a) RMSE derivatives (.)
(h) MNI.. (
10"
10 1
VPT
RMSE
10
10
2
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
SD of noise
SD of noise
(c) VPT (1)
(d) RMSE of energy (J)
Figure 10: Results with increasing standard deviation (SD) of Gaussian noise in training data.
Chua circuit
0.02
Ground Truth
HHD-GP
SPHHD-GP
0.01
Energy
0.00
-0.01
0.0
2.5
5.0
7.5
10.0
12.5
15.0
Time
Figure 11: Comparison of predicted energy of the Chua circuit.
K
Societal impact
Generally, our model can have a significant impact on a wide range of applications, such as phys- ical simulation and robotic control. However, it is imperative that meticulously evaluate the performance of prediction and uncertainty estimation when implementing it in societal contexts.
30

---

Ground TUB
Dourd Tul
StATPORE FURE
StatPoRE
HM5-0P Gread T2m
STIND-GP
-WW SHIP-2
SeriFort
400
Start Paret
Ground Stat
Grand
Stat Net
DW-SF Tull
Bal op Ground
AM Greans Trzm
1004 S1at AAT
Star PoRt
Stat
640
-40
Greand Trom
400
Groadd Tan
Figure 12: Comparison of predicted trajectories the Chua circuit.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The main contributions are summarized in the last paragraph of Section |
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: High computational complexity is the main limitation of our proposed model, which has been discussed in Appendix !
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 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 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 few datasets or with 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, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or 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 used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
32

---

Justification:
They are provided in Section 4 Section S and Appendix D
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibiliyy
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The information needed to reproduce the experimental results provided in Section Appendix Appendix F and Appendix G
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
33

---

Answer: [No]
Justification:
The code will be released afterwards.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https guwdes/dde/eeee Cooeeeoodeesssseeeeee 7/nips.com htd////bb///uuddddd//eeeeeieisss ssiniiiiccyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: The training and test details are provided in Section Appendix Appendix and Appendix G
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification:
We have provided the standard errors in Section
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
34

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
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 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 range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
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 (c.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http://eeurii..ccccpuccceetcccee
Answer:
[Yes]
Justification:
The research conducted in the paper conform with the NeurIPS Code of Ethics.
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 deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. 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 no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact 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.
35

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is 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 being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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 best faith effort.
12. Lieenses for existing assets
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?
Question:
Answer: [Yes]
Justification: The code packages and datasets are properly cited in the paper.
Guidelines:
The answer NA means that the paper does not 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 URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.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, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of 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.
36

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
37


Reviews:
[{"review_id": "HsWS55Bljp", "reviewer": "Reviewer_DJX1", "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.", "questions": "1.\tIs 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?\n2.\tOn line 328 it is stated that mirror symmetry is incorporated into the div-free kernel. Where is this kernel given?\n3.\tAre there general symmetries that should be held or is this always going to be problem/application specific? \n4.\tWhen the kernels do not admit a closed form how sensitive are they to approximations? And does this significantly affect the computational cost?\n5.\tIs there any impact of using sparse GPs? Ie are the sparse models still guaranteed to be G-invariant? \n6.\tRelevant GP citations:\n    o\tEquivariant Learning of Stochastic Fields: Gaussian Processes and Steerable Conditional Neural Processes, Holderrieth et al\n    o\tVector-valued Gaussian Processes on Riemannian Manifolds via Gauge Equivariant Projected Kernels – Hutchinson et al", "limitations": "yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is a well written, clear paper, that improves on recent work with convincing experiments.", "weaknesses": "1.\tThe 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.\n2.\tThe 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer DJX1\n\nWe are very grateful for your valuable comments and acknowledgement of our main contributions. We greatly appreciate the time and effort you put into reviewing our paper. Below are our responses to your questions.\n\n---\n\n> **Question 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?\n\n**Response**: The the Helmholtz-Hodge decomposition involves breaking down a vector field into its divergence-free and curl-free parts, which may not necessarily preserve the symmetries of the original vector field. This is because the decomposition process is based on mathematical properties (curl and divergence) rather than on the symmetries of the field. However, the symmetry-preserving Helmholtz-Hodge decomposition always exists for physical systems due to the physical relevance of their divergence-free and curl-free components. These physical relevance aligns well with the symmetries of underlying systems. The Helmholtz-Hodge decomposition respects these physical relevance, which naturally leads to the preservation of symmetries.\n\n> **Question 2**: On line 328 it is stated that mirror symmetry is incorporated into the div-free kernel. Where is this kernel given?\n\n**Response**: Do you mean the div-free kernel mentioned in line 378? We apologize that we did not provide details about this div-free kernel in the manuscript. To make our GP model identifiable, we constructed the div-free kernel by incorporating mirror symmetry with respect to the square boundary $[0, 1] \\\\times [0, 1]$. The symmetry group is given by \n$$\nG = \\\\left \\\\{ \n\\\\begin{bmatrix} 1 & 0 & 0\\\\\\\\ 0& 1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix}, \n\\\\begin{bmatrix} 1 & 0 & 0\\\\\\\\ 0& -1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 2\\\\\\\\ 0& -1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 2\\\\\\\\ 0& 1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 2\\\\\\\\ 0& -1 & 2\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} 1 & 0 & 0\\\\\\\\ 0& -1 & 2\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 0\\\\\\\\ 0& -1 & 2\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 0\\\\\\\\ 0& 1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix},\n\\\\begin{bmatrix} -1 & 0 & 0\\\\\\\\ 0& -1 & 0\\\\\\\\ 0 & 0 & 1 \\\\end{bmatrix}\n\\\\right \\\\}， \n$$\nwhere the group elements are the mirror transformation matrices applied to homogeneous coordinates $(x, y, 1)$. Theorem 5.3 in our manuscript shows that a $G$-equivariant div-free vector field can be constructed from a vector potential with the same symmetry. So, we constructed the kernel of the vector potential by using Eq.17 to integrate a *squared exponential kernel* over the symmetry group above. This potential kernel was then substituted into Eq.18 to construct the div-free kernel. The partial derivatives involve in Eq.18 were calculated by automatic differentiation in PyTorch, so we did not need to derive the analytic expression of the div-free kernel. We will add these details in the future manuscript. Thank you very much for pointing out this issue.\n\n\n> **Question 3**: Are there general symmetries that should be held or is this always going to be problem/application specific?\n\n**Response**: There is no general symmetry that must be maintained. Although we have indicated in the appendix that the actual required symmetry is problem-specific, we will make this more explicit in the main paper based on your suggestion.\n\n> **Question 4**: When the kernels do not admit a closed form, how sensitive are they to approximations? And does this significantly affect the computational cost?\n\n**Response**: The kernel maintaining rotation symmetry is constructed by integrating over the continuous rotation group, and this kernel has no closed form (Eq.72). So we numerically approximated this kernel by sampling discrete rotations. To investigate the impact of this kernel approximation, we added experiments on learning with varying number of sampling rotations. The results are given in the table below, from which we can observe that the performance of the kernel approximation improves rapidly with increasing number of sampling rotations. Kernels with rotational symmetry were well approximated even with a small number of samples (e.g., n = 8), which ensures that approximating the kernel by numerical integration does not significantly improve the computational cost. (The computational complexity increases linearly with the number of discrete rotation samples.)\n\n| Number of rotation samples | $n=1$ | $n=4$  | $n=8$  | $n=36$ | $n=100$ | \n|-----|-----|-----|-----|-----|-----|\n| RMSE $(\\\\times 10^{-2})$ |  5.76  |   3.12   |  1.96 | 1.74 | 1.67\n\n> **Question 5**: Is there any impact of using sparse GPs? Ie are the sparse models still guaranteed to be G-invariant?\n\n**Response**: In our method, the symmetry constraints are imposed to GP models by disigning suitable kernels. The sparse GP model we used reduces the computational complexity by using a smaller set of inducing points to compactly represent the original data, which do not alter the kernel's structure and properties. Therefore, the symmetries can still be respected in the sparse GP model.\n\n> **Question 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\n\n**Response**: Thank you for providing these two works, which will be cited in our future manuscript.\n\n---\nThank you again for the constructive comments. We hope these explanations could address your concerns. Any further questions or suggestions would be greatly appreciated.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "5PyR9Qe2b6", "reviewer": "Reviewer_w6Wk", "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.", "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].\n- Please add a discussion on computational complexity. Specifically, numerical integration is used to consider symmetry constraints, but could this be a computational bottleneck?\n- 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?\n\n[R1] Yusuke Tanaka, Tomoharu Iwata, Naonori Ueda, Symplectic Spectrum Gaussian Processes: Learning Hamiltonians from Noisy and Sparse Data, NeurIPS, 2022.", "limitations": "The differences between HHD-GP and existing methods (such as [R1]) should be thoroughly discussed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The curl-free and divergence-free kernel considering symmetry constraints is novel.\n- The manuscript is generally easy to read.\n- The model is evaluated on multiple dynamical systems.", "weaknesses": "- There are missing citations for relevant existing methods.\n- There is no discussion on computational complexity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer w6Wk\n\nThank you very much for your constructive comments and for taking the time to review our manuscript. We've carefully read your comments and our response is as follows.\n\n---\n\n> **Question 1**: 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].\n\n**Response**: SSGP and our models (HHD-GP and SPHHD-GP) study the problem of learning dynamical systems from different perspectives and scopes, leading to essential differences in their construction of kernels. Specifically, we emphasize that our kernels are constructed from the perspective of satisfying certain differential invariants (either free of divergence or of curl), rather than from physically governing equations (e.g. the Hamiltonian equation used by SSGP). Therefore, HHD-GP has a wider scope of applicability, which is applicable to any physical system that can be described by a smooth vector field, but SSGP can only be used for dissipative Hamiltonian systems. For example, SSGP cannot be used to model the Chua circuit, a chaotic system used in our experiments, because this system cannot be described by the Hamiltonian equation. Hamiltonian vector fields are a specific type of divergence-free vector fields. From this perspective, HHD-GP can indeed be considered as a generalization of SSGP. Furthermore, although SSGP combines a Hamiltonian kernel with an additive dissipation kernel, it does not consider the non-identifibility problem prevalent in additive models, making the model struggle with learning system components from noisy data. In contrast, we solve the non-identifibility problem in HHD-GP by incorporating symmetries of underlying dynamics. We believe this is an original and valuable contribution to the field of machine learning for dynamical systems. Thank you for providing us with the opportunity to further elaborate on the contributions of our work, which we will include in the future manuscript as you suggest.\n\n> **Question 2**: Please add a discussion on computational complexity. Specifically, numerical integration is used to consider symmetry constraints, but could this be a computational bottleneck?\n\n**Response**: When the symmetry-preserving kernel admits no closed form, numerical integration is exploited to approximate the kernel and compute the covariance, whose computational complexity grows linearly with the number of discrete points. To investigate the impact of the number of discrete points on model performance, we performed experiments on learning with the kernel maintaining rotation symmetry approximated by varying amount of discrete rotations. The results are given in the table below, from which we can observe that the performance of the kernel approximation improves rapidly with increasing number of discrete rotations. Kernels with rotational symmetry were well approximated even with a small number of samples (e.g., n = 8), which ensures that the numerical integration could not be a computational bottleneck. The computational complexity of our method is dominated by the Cholesky decomposition to invert the covariance matrix, which grows cubically with the number of training data points and the input dimension.\n\n| Number of rotation samples | $n=1$ | $n=4$  | $n=8$  | $n=36$ | $n=100$ | \n|-----|-----|-----|-----|-----|-----|\n| RMSE $(\\times 10^{-2})$ |  5.76  |   3.12   |  1.96 | 1.74 | 1.67\n\n\n> **Question 3**: 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?\n\n**Response**: Figure 3 provides the experimental results for the RMSE of energy prediction as the number of training data increases, which demonstrate the non-identifiability problem suffered by HHD-GP. The covariance matrix equation ($K_* = K_{curl} + K_{div} + \\\\Sigma$, Eq.6 in the manuscript) shows that when the HHD-GP is used to learn the decomposition, the effects of different components can be treated as observation noise to each other. However, the non-uniqueness of HHD makes the HHD-GP model non-identifiable, making it struggle to distinguish between the different factors contributing to the data. Therefore, more uncertainty in model parameters is introduced with more data, leading to an increase in prediction error as the number of training data increases. This inability of the model parameters to converge to their true values as the amount of data increases indicates a lack of model consistency, which is a necessary condition for a model to be identifiable, as discussed from line 321 of the manuscript. We will add more details to explain the experimental results, following your suggestion.\n\n---\nThank you again for the constructive comments. We hope these explanations could address your concerns. Any further questions or suggestions would be greatly appreciated.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "F8dGA0Jww9", "reviewer": "Reviewer_X44d", "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.", "questions": "-\tPlease 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).\n-\tRelated 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": "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.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "-\tThe approach taken is novel, but not unnecessarily complex. The theoretical insight regarding the Helmholtz-Hodge decomposition is combined with a basic building block of GP-based modelling (additive combination of GPs) in order to yield a technique which consistently outperforms all the relevant baselines presented by the authors.\n-\tThe quality of the work carried out in general is very good; the experiments are thorough and well documented overall. The Hamiltonian dynamics experiments are effectively a standard benchmark in this area, the chaotic system presents a further challenge for the model, and the ocean current experiment is a real-world example of how this approach may be useful in a practical setting. The performance of the SPHHD-GP across all of these experiments is impressive.\n-\tThe clarity of presentation is also very good, and I believe that the work would be accessible both to researchers from a dynamical systems background without in-depth knowledge of GPs, and vice versa. The narrative of the paper is well crafted, from problem setting to initial formulation of the model, onto discussion of limitations regarding identifiability, and finally the implementation of symmetry-based constraints which address these limitations.", "weaknesses": "-\tFrom 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.\n-\tSeveral 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer X44d\n\nThank you very much for your thoughtful review and constructive comments. We are very pleased that you recognize the importance and contribution of our work. We have carefully gone through your comments and suggestions, and we believe addressing these points in the manuscript would indeed make the paper better.\n\n---\n\n> **Question 1**: 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).\n\n**Response**: Following your suggestion, we will add some brief details about the implementation of the sparse GP into the appendix. The key idea of this sparse GP framework is to introduce pseudo-inputs, which form a smaller set of data serving as a compact representation of the original data. These pseudo-inputs are not necessarily a subset of the training data but are optimized to capture the essential structure of the dataset. In our experiments in Section 6.3, 200 training data points are randomly selected as the initial pseudo-inputs, whose locations are learned jointly with kernel parameters through the ADAM optimizer. These pseudo-inputs greatly reduce the computational complexity while ensuring that the sparse GP model maintains a high prediction accuracy.\n\n> **Question 2**: 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)?\n\n**Response**: Considering the computational feasibility, the experiments in Sections 6.1 and 6.2 were performed using exact inference. The sparse GP framework can achieve similar prediction accuracy to the exact inference as long as sufficient pseudo-inputs are employed.\n\n---\nThank you again for the constructive comments. We hope these explanations could address your concerns. Any further questions or suggestions would be greatly appreciated.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "J7WkxBlTmC", "reviewer": "Reviewer_KheP", "summary": "This paper proposes a novel dimensionality reduction method.\nThe method relies on the Helmholz-Hodge decomposition to identify the dynamical system through a decomposition into a curl-free and a divergence-free part.\nFurthermore, the method introduces a way to incorporate priors that constrain the identified vector fields based on symmetry constraints.\nThis method is applied to various toy models (mass-spring system, pendulum, Chua circuit) and to real world data (ocean current fields).", "questions": "What is the relevance of $\\mathcal{G}$-equivariance and uniqueness/identifiability?\nAppendix H: Uniqueness and Symmetries\n- How do you know what priors to choose? Isn't an agnostic method better?\n- How do you know which symmetries to enforce?\n\nCould it be more fruitful to restrict yourself to a compact manifold with boundary so that the decomposition is unique?\n\n\n\n# Implementation details\n- How were the kernel parameters and noise variances initialized? Does training depend on initialization?\n(Here you mention that these are hyperparameters, but typically one refers to non trained parameters as hyperparamters.)\n-What were the ADAM optimizer paramters used for training?\n\n\n# Choice of evaluation metrics \n-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.\nWhy not prediction of the future state/trajectory? Why prediction only measured through VPT?", "limitations": "The authors should comment on limitations and societal impact.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper addresses an important problem in dynamical systems modeling.\nThis paper is the first to introduce the use of symmetries to construct unique GP models of dynamics that are decomposed into div- and curl-free vector fields.\nThis method is shown to be consistent in terms of convergence to the true parameter values as the amount of data increases.\nThe paper discusses the computational complexity of the proposed method.\nThe paper demonstrates that the presented framework is outperforms baselines on various toy models\nFurthermore, the method is shown to be able to handle noisy data.\nFinally, the method is applied to a real world dataset, which indicates practical significance.", "weaknesses": "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.\n\nThe paper mentions that the accurate prediction of energy demonstrates the interpretability of learned div-free features. \nHowever, this is not well-explained. The paper claims that the curl- and div-free components are physically meaningful, but this should be better justified.\n\nThe paper is missing comparison to other existing approaches for modeling physical systems. See for example: \nGuo, Q., Mandal, M. K., \\& Li, M. Y. (2005). Efficient Hodge–Helmholtz decomposition of motion fields. Pattern Recognition Letters, 26(4), 493-501.\n\nCompare to other dynamical systems reconstruction methods:\n- Hess, F., Monfared, Z., Brenner, M., \\& Durstewitz, D. (2023). Generalized teacher forcing for learning chaotic dynamics. arXiv preprint arXiv:2306.04406.\n\n-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.\n\nFinally, there is no code provided to show the method's working and not all hyperparameters are described (see also Questions).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Response to Reviewer KheP\n\nWe are very grateful for your valuable comments and acknowledgement of our main contributions. We sincerely appreciate your time and effort in reviewing our paper. Here is our response to your comments.\n\n---\n\n> **Question 1**: 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?\n\n**Response**: \n\nThe current work assumes that the prior of symmetries is directly available, and focuses on developing the method of incorperating symmetries into a GP model. This assumption is reasonable for many dynamical systems, but we totally agree with your concern about the availability of symmetry. Although physical systems always adhere to symmetries, the symmetries of many complex phenomena are often considered to be partially known or unknown. Now that our model has been demonstrated to perform well when system symmetry is available, we look forward to developing a method for symmetry-agnostic scenarios, i.e. enhancing our model with the capabilities of automatedly learning symmetries. Several recent studies have explored the feasibility of learning hidden symmetries from data [1, 2]. Since our model is developed in a probabilistic framework, another promising direction is to incorporate approximate symmetry instead. Studying the relationship between approximate symmetry and model identifiability is definitely an interesting and promising but unexplored research problem.\n\nWe also agree with you that enforcing boundary condition is another effective way to make the decomposition unique. Our method of enforcing symmetry can also be used to impose the boundary condition, which is demonstrated by our experiments on learning ocean current fields (Section 6.3 of the manuscript). In this experiment, mirror symmetry was incorporated into the div-free kernel to enforce the divergence-free vector field to be parallel to the domain boundary, which is a sufficient condition for the decomposition to be unique [3]. The experimental results in Section 6.3 of the manuscript demonstrate that enforcing the boundary condition allows our model to provide more realistic ocean current predictions and divergence identification than the baselines.\n\n> **Question 2**: How were the kernel parameters and noise variances initialized? Does training depend on initialization? What were the ADAM optimizer paramters used for training?\n\n**Response**: For the experiments on simulated systems (damped mass-spring system, damped pendulum, and Chua circuit), we set the noise level to 0.05 (lines 275 and 355 in the main text). To explore the impact of noise level on model performance, in Appendix J.5 we provide experimental results with increasing noise in the training data (0.01, 0.05, 0.1, 0.2). For the experiments on a realistic system (ocean current field), we did not manually add noise to the training data since we used a realistic dataset. The kernel parameters are initialized randomly from a range of [0.1, 10]. In Table 1 of the manuscript we reported experimental results averaged over 10 independent experiments performed by resampling the kernel initial parameters, from which we observed no dependence of training on the initial parameters. As discussed in line 802, the GP-based models were trained by the ADAM optimizer, with a learning rate of 0.01 for 3000 gradient steps, and other parameters are the default values ​​in pytorch.\n\n> **Question 3**: 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?\n\n**Response**: In our experiments, we evaluate the performance of the models in predicting both vector fields and state trajectories, with the former measured by the root mean squared error (RMSE) and the latter measured by the valid prediction time (VPT). The reason why we chose VPT is that many recent studies in learning dynamical systems [4,5,6] suggested that the RMSEs of state trajectories over long time horizons can be misleading indicators. As an example in the pendulum dataset, a trajectory remaining stationary at its initial position may have lower RMSE compared to a trajectory recovering the oscillatory behavior correctly but having a slight shift in angular velocity. Therefore, VPT has become a popular metric for measuring a model's ability to predict state trajectories.\n\n---\nThank you again for the constructive comments. We hope these explanations could address your concerns. Any further questions or suggestions would be greatly appreciated.\n\n---\n## References\n\n[1] Liu, Ziming, and Max Tegmark. \"Machine learning hidden symmetries.\" Physical Review Letters 128.18 (2022): 180201.\n\n[2] Desai, Krish, Benjamin Nachman, and Jesse Thaler. \"Symmetry discovery with deep learning.\" Physical Review D 105.9 (2022): 096031.\n\n[3] Bhatia, Harsh, et al. \"The Helmholtz-Hodge decomposition—a survey.\" IEEE Transactions on visualization and computer graphics 19.8 (2012): 1386-1404.\n\n[4] Matsubara, Takashi, and Takaharu Yaguchi. \"Finde: Neural differential equations for finding and preserving invariant quantities.\" arXiv preprint arXiv:2210.00272 (2022).\n\n[5] Jin, Pengzhan, et al. \"SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems.\" Neural Networks 132 (2020): 166-179.\n\n[6] Vlachas, Pantelis-Rafail, et al. \"Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.\" Neural Networks 126 (2020): 191-217.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: 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.
Decision: Accept (poster)