Symmetry-Informed Governing Equation Discovery
Authors: Jianke Yang
Wang Rao
Nima Dehmamy
Robin Walters
Rose Yu
Venue: NeurIPS 2024 poster (2024)
Keywords: Equation Discovery
Dynamical Systems
Symmetry
Equivariant Neural Networks
Symbolic Regression
Lie Group
Ordinary Differential Equation
Lie Point Symmetry
Primary area: machine_learning_for_physical_sciences
TL;DR: A new approach to improve equation discovery of dynamical systems with symmetry priors

Abstract: Despite the advancements in learning governing differential equations from observations of dynamical systems, data-driven methods are often unaware of fundamental physical laws, such as frame invariance. As a result, these algorithms may search an unnecessarily large space and discover less accurate or overly complex equations. In this paper, we propose to leverage symmetry in automated equation discovery to compress the equation search space and improve the accuracy and simplicity of the learned equations. Specifically, we derive equivariance constraints from the time-independent symmetries of ODEs. Depending on the types of symmetries, we develop a pipeline for incorporating symmetry constraints into various equation discovery algorithms, including sparse regression and genetic programming. In experiments across diverse dynamical systems, our approach demonstrates better robustness against noise and recovers governing equations with significantly higher probability than baselines without symmetry.

Symmetry-Informed Governing Equation Discovery
Jianke Yang UCSD
Wang Rao" Tsinghua University
Nima Dehmamy IBM Research
Robin Walters Northeastern University
Rose Yu UCSD
Abstract
Despite the advancements in learning governing differential equations from observa- tions of dynamical systems, data-driven methods are often unaware of fundamental physical laws, such as frame invariance. As result, these algorithms may search an unnecessarily large space and discover less accurate or overly complex equations. In this paper, we propose to leverage symmetry in automated equation discovery to compress the equation search space and improve the accuracy and simplicity the learned equations. Specifically, we derive equivariance constraints from the indeeindeeendent symmetries of ODEs. Depending on the types of symmetries, develop pipeline for incorporating symmetry constraints into various equation discovery algorithms, including sparse regression and genetic programming. experiments across diverse dynamical systems, our approach demonstrates bet- ter robustness against noise and recovers governing equations with significantly higher probability than baselines without symmetry. Our codebase is available at https::/githu..com/RRe
Ltwtaa/o/oooooo-S-Laaaaaeeeeeeee STLL-Labyymeety--oe----ee
1 Introduction
Discovering governing equations that describe and predict observational data from nature is primary goal of science. Compared to the black-box function approximators common in machine learning, symbolic equations reveal deeper understanding of underlying physical processes. In this paper, we consider a first-order dynamical system governed by an autonomous ordinary differential equation (ODE):
i(t) = hxx(t))
(1)
where E T denotes time, xt) E X is the system state at time t, X C IRd is the phase space of the ODE, (t)) the time derivative, and h Rd the vector field describing the evolution dynamics of the system. Our goal is to discover the function h in a symbolic form from the observed trajectories of the system.
Traditionally, discovering governing equations has been difficult task for human experts due to the complexity of the underlying physical systems and the large search space for equations. Recently, many new algorithms have been proposed for automatically discovering equations from data. When applied an ODE, these methods first estimate the label, i.e. the time derivative, by numerically differentiating the trajectory. A variety of techniques including sparse function basis regression (Brunton et al. 2016 Champion et al., 2019) and genetic programming (Cranmer et al., 2019) are then used to fit the function between the measurement and the estimated derivative. Another class of methods considers the differential equation's variational form eliminate the need for derivative estimation and improve the robustness against noise (Messenger & Bortz 2021
Qian et al. 2022).
An important aspect of equation discovery is the principle of parsimony (Sober) 1981 the scientific principle that the most acceptable explanation of phenomenon is the simplest. It promotes the scientific goal of identifying a model that offers robust predictive capabilities while maintaining simplicity. Existing methods that fail to consider these principles may end up searching in an
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

unnecessarily large space of equations and discover equations that are overly complex or do not conform fundamental physical laws.
In this work, we advocate for using the principle of symmetry to guide the equation discovery process. Symmetry has been widely used in machine learning design equivariant networks (Bronstein et al.)
2021 Finzi et al. [222 Wang et al. (022), augment training data (Benton et al. (020), or regularize
the model AAkoonddSSddeh et al. 2023] Otto al. 2023), leading improved prediction accuracy,
sample efficiency, and generalization.
Specifically, derive the equivariance constraints from the symmetries of ODEs and solve the constraints explicitly to compress the search space of equations. Alternatiely, we can also use symmetry as regularization loss term to guide the equation discovery process and improve its robustness to measurement noise. In practice, a dynamical system may possess symmetries that are unknown or difficult to describe in simple terms. To account for these situations, also incorporate method for learning symmetries with nonlinear group actions (Yang et al.) (203a)). We then establish a pipeline capable of learning unknown symmetries from data and subsequently discovering governing equations that conform to the discovered symmetry. To summarize, our main contributions include:
We establish a holistic pipeline to use Lie point symmetries of ODEs to improve the accuracy and robustness of equation discovery algorithms.
We derive the theoretical criterion for symmetry of Lime-independent ODEs in terms of equivariance of the associated flow map.
From this criterion, we solve the linear symmetry constraint explicitly for compressing the equation search space in sparse regression and promoting parsimony in the learned equations.
For general symmetry, we promote symmetry by symmetry regularization term when the symmetry constraint cannot be explicitly solved.
In experiments across many dynamical systems with substantial noise, our symmetry- informed approach achieves higher success rates in recovering the governing equations and better robustness against noise.
2 Related Works
2.1 Governing Equation Discovery
There are numerous methods for discovering governing equations from data. Genetic programming has been successful in searching exponentially large spaces for combinations of mathematical operations and functions (Gaucel et al.) 2014). Recent works have applied genetic programming to distill graph neural networks into symbolic expressions (Cranmer et al. 2019 2020). Udrescu & Tegmark| (2020) introduce physical inductive biases to expedite the search. The main limitation of genetic programming arises from computational complexity due to the combinatorially large search space. Also, genetic programming can suffer from noisy labels (Agapitos et al. 2012). which are common in discovering ODEs due to measurement noise.
Another branch of research involves sparse regression. Originating from SINDy (Brunton al. 2016). these methods assume that the equation can be written linear combination of some predefined functions and optimize the coefficients with sparsity regularization. Subsequently, [Chmmpion et al.) (2019) combines the sparse regression technique with an autoencoder network to simultaneously discover coordinate systems and governing equations. Weak SINDy (Messenger & Bortz, 202T) uses an alternative optimization objective based on the variational form of the ODE, which eliminates the need for estimating time derivatives and improves the robustness against measurement noise. Other works improve upon SINDy by incorporating physical priors, such dimensional analysis (Xie et al. 2022 Bakarji et al. 2022). physical structure embedding (Lee et al. 2022). and symmetries (Loiseau
& Brunton 2018 Guan et al. 2021.
2.2 Symmetry Prior in Machine Learning
Symmetry plays erucial role in machine learning. There has been a vast body of work on equivariant neural networks (Cohen & Welling 2016; Weiler & Cesa. 2019: Finzi et al. 2020 2021 Bronstein

---

et a.. 2021 Wang et al. 2021). These network architectures enforce various symmetries in different data types. Symmetry can also be exploited through data augmentation (Benton et al. 2020) and canonicalization (Kaba et al. 2023). There are already some works that exploit symmetries for learn- ing dynamical systems governed by differential equations through data augmentation (Brandstetter et al., 2 symmetry loss (Huh et al. 2020 Akhound-Sadegl. et al.) 2023) and contrastive learning (Mialon et al. 2023). However, works using symmetry in recovering underlying equations from data are scarce, and existing examples only consider symmetries specific to a particular system, e.g. reflections and permutations in proper orthogonal decomposition (POD) coefficients (Guan et al.
202Ib, rotations and translations in space (Ridderbusch al. 2021 Baddoo et al. 2021 etc. In this work, we provide solutions for dealing with general matrix Lie group symmetries.
While symmetry has proved an important inductive bias, for the ODEs considered in this work, the underlying symmetries are often not known priori. Some recent works have aimed to automatically discover symmetries from data (Liu & Tegmark [022)] Yang et al. 2023b a Otto et al. [023] van der Ouderaa et al. 2024). When there no available prior knowledge about symmetry. We use the adversarial framework in Yang al. (2023a) to learn unknown symmetries from data and discover the equations subsequently using the learned symmetry.
3 Symmetry of ODEs
Our goal is to discover the governing equations of a dynamical system (T) in symbolic form. Our J/Y
dataset consists of observed trajectories of the system, { where at X Rd denotes the observation at timestep t in the ith trajectory. We assume the observations are collected at regular time interval At. We omit the trajectory index and consider one trajectory. To:T for simplicity. To leverage symmetry for better equation discovery, we first define symmetries in ODEs and provide the theoretical foundation for our methodology. In particular, we show that -nneenndeenddet symmetry of the ODE is equivalent to the equivariance of its flow map. This equivariance property can then be conveniently exploited to enforce symmetry constraints in the equation discovery process.
3.1 Lie Point Symmetry and Flow Map Equivariance
The Lie point symmetry of an ODE is transformation that maps of its solutions another. The symmetry transformations form the symmetry group of the ODE. In this paper, we focus on symmetry transformations that solely on the phase space X æ. without changing the independent variable I. We refer to these as time-independent symmetries, formally defined as follows.
Definition 3.1. A iimeinndeenndent symmetry group of the ODE IT a group G acting on X such that whenever x: is a solution of ( x g ax()) is also solution of the system.
We refer the readers to Appendix A for background on Lie groups, group actions, and their applications on differential equations. We consider time-independent symmetries in Definition 3.1 and demonstrate how to use these symmetries derive equivariance constraints on equations.
Given fixed time т, the flow map associated with the ODE ( is denoted fr : X, which maps starting point in trajectory æ(to) to endpoint ætto + T) after moving along the vector field h for time T. Formally, f=o)) = x(T). where a+()) is the solution of the initial value problem
a (t) hxx,, x x0)) = zo. t [0.
(2)
Proposition 3.2. Let G be a group that acts on the phase space X of the ODE Q G is symmetry group of the ODE ① in terms of Definition 3.1 if and only if for any T E T, the flow map fr is
equivariant to the G-action on X.
As our dataset consists of trajectories measured at discrete timesteps, consider the discretized flow map 2++1 f(x)). where f == fa. depends the sampling step size At of the data. From Proposition 3.2 this function of moving one step forward in the trajectory is G-equivariant, i.e.
Ut+1 f(x)) = g It+1 x()
(3)
3

---

3.2 The Infinitesimal Formulae
We can use the equivariance condition 3 constrain the equation learning problem. Consider the Lie group of symmetry transformations G and the associated Lie algebra g. We have the following equality constraints.
Theorem 3.3. Let G be a time-independeni symmetry group of the ODE ( Let v be an element in the Lie algebra g of the group G. Consider the flow map f defined 22for a fixed time interval. Denote J and Ig as the Jacobian of fx) and g a w.E.t a. For all x E X and g = exp(ev) E G, the following equalities hold:
f(g x) 9 f(x)) 0;
(4) (6)
))( - v(f())) 0; ( J.(x))(()) 0. (7)
(5)
(x))) h0g x) O;
To understand Theorem 3.3 follows directly from Proposition 3.2 and the other equations are the infinitesimal equivalent of this equivariance condition. For example, the Jacobian-vector product in (5) reveals how much f changes when the input infinitesimally transformed by v, and v v (a)) indicates the amount of change caused by the infinitesimal action on the output of f. F. The difference between these two terms vanishes if the function equivariant. On the other hand, as both h and can be viewed as vector fields on X, their roles are interchangeable in the equivariance formula. Therefore, we can also express the equivariance in terms of the flow h and a finite group element g exp(ev). Finally, we also consider the fully infinitesimal formula (7). which indicates that h and v commute as vector fields. The proof of Theorem 3.3 may be found in Appendix [.I
Theorem provides the theoretic basis for our symmetry-basee based equation discovery methods in Section We can enforce the equality constraints in Theorem 3.3 by solving the constraints explicitly (Section or adding them as regularization terms (Section
4 Symmetry-Informed Governing Equation Discovery
Next, we discuss how leverage the knowledge of symmetry in equation discovery. Depending on the nature of different symmetries, we develop a holistic pipeline for incorporating symmetries into various equation discovery algorithms, shown in Figure following subsections.
We introduce the components in the
4.1 Equiv-c: Solving the Linear Symmetry Constraint
We start from the case of linear symmetries. Many ODEs exhibit symmetries with linear actions on the state variable æ. For example, in classical mechanics, the motion of a particle under a central force is invariant under rotations about the center. For such linear symmetries, if the governing equations can be written as linear combinations basis functions as in SINDy sparse regression (Brunton et al.. 2016), the infinitesimal formula 77) becomes linear constraint. The constraint can be explicitly solved to construct an equivariant model for h.
In particular, suppose the governing equations can be written in terms of O(x) RP as:
h(x)= (
(8)
where W E Rdxp is the learnable coefficient matrix. In SINDy, the function library ((x)) is pre- defined. For example, @ can be the set of all polynomials up to second order for a system, i.c. ((11, x2) = [L, æ1, x2, 2f 2132, x5
2 2-dimensional
Our goal is to solve for the parameter W that makes 7 hold for any æ. To this end, we define the symbolic map Me (IR' + R) - IR" which maps a function to its coordinate in the function space spanned by Ô. Specifically, Mf(
for fz( c Oz)), c R.. For a multivariate function f: Rdd + F., we compute Me for each of its components and stack them into Me()) E Raxp, A concrete example is provided in Appendix B.2 We have the following proposition.
Proposition 4.1. Let G be a time innepenndeennentt symmetry group of the ODE ē W@(x) with linear actions. Let
g be Lie algebra element with action X by vi(ar) T,yær T.p E Rdxd. Then, L. = W Me Me(Je (.). where Je denotes the Jacobian of O.
Note that v( v((())) means evaluating the vector field at f(x) E X. is different from v (z)) which is the Lie derivative of at the point

---

Data .07
Symmetry
Equation learning j h(()
Knowledgge about symmetry?
No
Discover symmetry
a
Discover Smmmerr
Symmeriy regullrization
[[[v. A]]l
Linear symmetry?
No
Symmetry regularization
by sunnn
b
Yes
No
Q: Equivariant basis
Learnable parameters
Linear equation form?
Latent equation?
symmetry
EquisSINDY
@
h(()) = w arg min, || 1(x0)12
EquivSINDy +AE
Enc.
EquivSINDY
C
C
G
Figure 1: Pipeline for incorporating symmetries into equation discovery via solving linear symmetry constraint (Section 4.1), regularization (Section 4.2) and symmetry discovery (Section 4.3). Given the trajectory data from the dynamical system, we first identify its symmetry based on prior knowledge symmetry discovery techniques. We then enforce the symmetry by solving set of constraints when possible and otherwise promote the symmetry through regularization.
Proof. Substituting 8 into 77 we obtain J. (O))) = ()))())((x), Vo E X. As acts on X linearly, i.e. v(x)= L, the constraint equivalent to L.O W (e))
We then apply Me to both sides the equation. Obviously, Me(O) = Ip. Since is true for all æ, have L.W. W = W (JMe(Ja()Lo.))
To calculate Me, we need to ensure that (e) is still within the span of the functions in the given (x:). This is true when ((x)) is the of all polynomials up certain degree, which prove in Appendix [B.2]
Proposition 4.2. The components of Jo L,, E RR be written linear combinations the terms in (x)) if (x:) is the set of all polynomials up to degree E z
Once we have the infinitesimal action Ty and appropriate function library θ, Me can be computed in purely symbolic way using sympy (Meurer al., 2017). Denoting Me L. = Me (., and vectorizing the matrix W, we can rewrite the equation Proposition 4.1 as the following linear constraint on W:
(- LL,)voc =
where ① is the Kronecker sum: A@B = A & I + I & B. We concatenate the constraints for all representations I.,, of the Lie algebra basis {u;}g_) into a single matrix C and solve the constraint using the singular value decomposition:
M Mb.L., BLm
C vec(W)=
vec(W)
MM. BLv,
rec(W) = 0.
(10)
0
The coefficient matrix resides in the r-dimensional nullspace of C and can parametrized as vec(W) = Q3, 3 E IR. This equivariant parametrization significantly reduces the number of free parameters from x p to r, leading to parsimonious equations and easier training compressed parameter space. We name this approach EquivSINDy-c where stands for constraint.

---

The above procedure uses linear symmetries to inform sparse regression. However, when the data is nigh-dimensional, such as videos, linear symmetries can become too restrictive. For these scenarios, SINDy autoencoder (Champion et al., 2019) maps the data to a latent space to discover the coordinate system and the governing equation under that system. We show in the experiment section that our model can also be used in this case by enforcing the latent equation to be equivariant. With the inductive bias of symmetry, the resulting coordinate and equation can achieve better long-term prediction accuracy.
4.2 Equiv-r: Regularization for General Symmetry
The above procedure, Equiv applies to linear symmetry with a proper choice of function library é. Solving the equality constraint in Proposition [.1 for arbitrary group action is generally more challenging. Technically, when the infinitesimal action is expressed in closed form, we can still apply Me to both sides of the equation and obtain Ma(r(WO(()) = W (Me (J)())u())) However, the function library needs to be carefully chosen based on the specific action v to ensure that Me is evaluated on functions within the span of (az). Moreover, as we will discuss in Section 4.3 we may rely on symmetry discovery methods to learn unknown symmetries parametrized by neural networks. Without a closed-form symmetry, calculating Me is intractable.
We propose alternative approach that is universally applicable any symmetry and equation discovery algorithm. We use the formula from Theorem 3.3 as a regularization term to promote symmetry in the learned equation. This symmetry loss is added the equation loss from the base equation discovery algorithm, c.g. L2 error between the estimated time derivative and the prediction from the learned equation.
We consider the following relative loss based on infinitesimal group action e and a finite-time flow map fr of the equation in
|Jg, (z))(()) - u(f
Cuymu
Eg
Σ (E
(11)
[
where Bgg) denotes the basis of the Lie algebra, f- is obtained by solving the initial value problem (2) and Jy, denotes the Jacobian of fr.
We use a relative loss because the scale of both terms in the numerator and their difference subject to the specific fr. In the extreme case when h O, both terms are zero because (x) is constant. As we are optimizing h under this objective, we compute the relative error to eliminate the influence of the level of variations in the function itself on the symmetry error.
We can also introduce other forms of symmetry regularizations based on the other equations in Theorem 3.3 Empirically, these regularization terms perform similarly, but their implementations have subtle differences. For instance, some loss terms require integrating the learned equations to get the flow map f- and thus incur additional computational cost; some involve computing higher-order derivatives for infinitesimal symmetries, which amplifies the numerical error when the symmetry is not exactly accurate. We list these different options of symmetry regularization and compare them in more detail in Appendix C.2
The use of symmetry regularization encourages equation discovery models to achieve lower symmetry error, instead of explicitly constraining the parameter space as in Section 4.1 However, is general approach that can be applied to various equation discovery algorithms, e.g. sparse regression and genetic programming. Similar to the works that utilize the variational form of ODEs (Messenger & Bortz 2021: Qian et al. (2022), , our symmetry regularization does not require estimating the time derivative, so it is more robust to noise. We name this approach as EquivSINDy-r (or EquivGP-r for genetic programming- based discovery algorithms), standing for regularization.
4.3 Equation Discovery with Unknown Symmetry
The knowledge of ODE symmetries often not accessible when we do not know the equations. In this case, we can use symmetry discovery methods (Dehmamy et al.) 2021 Liu & Tegmark 2022] Yang al. 2023b a Otto et al., (023) to first learn the symmetry from data and then use the discovered symmetries to improve equation discovery.

---

In our experiments, we use LaLiGAN Yang et al. 2023a) to learn nonlinear actions of a Lie group G I GL(K; as exp(e0)æ = (v exp(eL,,) ) o (z), where the networks v and define an autoencoder and Lg E R** is the matrix representation of E g. Our dataset consists of ODE trajectories sampled at a fixed rate At, each given by 20.T. We can learn the equivariance of 3t 2c+1 by feeding the imput-output pairs (0;, into LaLiGAN. We extract the infinitesimal action as
d (a)) = de
(vb o exp(ee.. ()()) ==0
(12)
where z == =x)) is the latent mapping via encoder network, and the Jacobian- vecttor product can be obtained by automatic differentiation. This infinitesimal action is then used to compute the regularization (T].
Alternatively, we can also discover equations in the latent space of LaLiGAN by solving the equiv- ariance constraints in Section 4.1 as the latent dynamics 2 + 2(+1 has linear symmetry L Discovering governing equations in the latent space can be helpful when the dataset consists of high-dimensional observations instead of low-dimensional state variables.
We should note that our method has a completely different goal from LaLiGAN. Our method aims inductive bias. LaLiGAN aims to discover unknown
to discover equations using symmetry symmetry. Only when we do not know the symmetry a priori, we use LaLiGAN as a tool to discover the symmetry first and use the discovered symmetry regularize equation learning.
Experiments
The experiment section is organized as follows. In Section 5.1 we will demonstrate how solve the linear symmetry constraints (Section[4.) on some synthetic equations. In Section [5.2] we will apply the same equivariant model discover equations in symmetric latent space for some high- dimensional observations. In Section [.3] we consider two well-studied benchmark problems for dynamical system inference under significant noise levels, where use the regularization technique (Section .2 based on learned symmetries.
Data generation. For each ODE system, we sample N random initial conditions æn from a uniform distribution on a specified subset of X C Rd. Starting from each initial condition, we integrate the ODE using the 4th-order Runge-Kutte (RK4) method and sample with a regular step size Al to get the discrete trajectory To.T- We use different N, T and At for each system, as specified in Appendix
Then, unless otherwise specified, we add white noise to each dimension of the state æ at a noise level d-- The scale of the noise depends on the variance in the data within each state dimension my; o; RR (z;). We will report the noise level dn for each system. After adding the noise, we apply a Gaussian process-baat based smoothing procedure on the noisy data, similar to Qian et al. (2022).
Evaluation metrics. We run each algorithm multiple times and report the success probability of discovering the correct function form. Specifically, assume the true equation is expanded as E (xx),, where fy: is nonzero constant parameter and f; is function of æ. Also, expand the discovered equation as E0 f((r), where 0;. * 0. The discovery is considered successful all the terms match, i.e. {ff} = {.f:). The probability is computed as the proportion of successful runs. As each dynamical system consists of multiple equations (one for each dimension), the success probability of discovering each individual equation (e.g. "Eq. 1". "Eq. 2" columns in Table T and the joint probability of successfully discovering all the equations (c.g. "All" column in Table T. We argue that this is the most important metric for evaluating the performance of an equation discovery algorithm, since an accurate equation form reveals the key variables and their interactions, and thus the underlying structure and relationships governing the system.
evaluate
The following metrics are also considered:
RMSE for parameter estimation Σ 0 @(h) |IZK, where K is the number of runs, 0 (0,, fa...) is the constant parameters in true equation terms defined above, and
7

---

( the estimated parameters of the corresponding terms from the kth run. If fi is not in the learned equation, d,. O. We evaluate this metric on (1) all runs and (2) successful runs, referred to as RMSE (all) and RMSE (successful). Also, for each system, this metric is computed for each individual equation and all equations.
Long-term prediction error. We integrate the learned ODE from an initial condition sampled from the test predict the future trajectory. We evaluate the error between the predicted and the true trajectory at chosen timesteps.
Algorithms and Baselines. For baseline comparisons, we include SINDy sparse regression (Brun- ton et al. 20T6) and Genetic Programming (GP) with PySR package (Cranmer 2023), and their respective variants using the variational formulation, i.e. Weak SINDy (WSINDy) (Messenger & Bortz 202I) and D-CODE (Qian et al., 2022). Based on different experiment setups and the guide- lines in Figure D we apply our methods referred to EquivSINDy /
EquivGP, where the suffix stands for solving the linear symmetry constraint, and r stands for using the symmetry regularization.
5.1 Equivariant SINDy for Linear Symmetries
Damped Oscillator: Qz, Q2
In this section, we consider two dynamical systems which possess some linear symmetries. These sym- metries are commonly observed across many systems and their presence can be easily detected. Therefore, we assume we know these symmetries and focus on symmetry-informer equation discovery by solving the symmetry constraints.
Growth: Q1.2,
Aj4
Ai
Damped Oscillator. This is two-dimensional sys- tem with rotation symmetry, characterized by the in-
Figure 2: Solutions for equivariant con- straints of (13 [4). In (13), the 2D param- eter space is spanned by Q1 and Q2. In (14). the 3D parameter space is spanned by
finitesimal generator e x201 01 (i.e.
-1 0]
in matrix terms). It is governed by 13
Q1,2,3, each marked with different color.
d 0.12 x2 db = 21 0.1æ²
dn = 0.33m + 0.1
(13)
(14)
42
Growth. This system exhibits scaling symmetry under (21, + (a²x,, axz). The correspond- ing infinitesimal generator v = 2x, 0 + x2222 is governed by [4).
Figure visualizes the parameter spaces for these two systems under the equivari- ance constraint. The original parameter build the func- tion library è with to second-order polynomials. For (13), the equivariance constraint reduces the parameter space two dimensions. The equivariant basis Q Qei and Q2 Qe2 are visualized. For (14), the parameter space is reduced to three dimensions. Each equivariant ba- sis component only affects one equation term, visualized in three distinct colors. We can observe that the equivariant mod- els have much fewer free parameters than non-equivariant ones.
Table 1: Success probability of equation discovery on the damped oscillator (3) at noise level o = 20% and the growth system [4) at 5%, computed from 100 runs for each algorithm. See Appendix C.T for full results.
space is W _22
Success prob.
System
Method
Eq.
Eq.
All
GP
0.00
0.70
0.00
D-CODE
0.00
0.00
0.00
Oscillator
13
SINDY
0.35
0.38
0.15
WSINDY
0.06
0.07
0.00
EquivSINDy-c
0.93
0.97
0.90
GP
0.00
1.00
0.00
D-CODE
0,00
0.65
0.00
Growth 14
SINDY
0.26
0.13
0.03
WSINDY 0.00
0.00
0.00
EquivSINDy-c
1.00
1.00
1.00
Table displays the success probability of different methods. With the reduced weight space, our equivariant model almost always recovers the correct equation forms. Full results, including the parameter estimation error and the prediction error with discovered equations, as well as the experiment on a higher-dimensional system, are deferred to Appendix C.T Overall, these results fully demonstrate the advantage of applying symmetry to compress the parameter space.
8

---

5.2 Equivariant SINDy in Latent Space
We also show an example of discovering equa- tions in the latent space with an equivariant model. We consider À w reaction-diffusion system visualized in Figure 3 The system is sampled at a 100 100 spatial grid, yielding observations
Prediction Ertor on Beactian-Difluiion
800 (og sca90
R 10000
See Appendix D.I for more details.
We apply LaLiGAN YYang et al., 2023a) dis- cover a latent space R where the dynamics 2t + Zt+1 is equivariant to linear group ac- tion. Meanwhile, we use equivariant SINDy with the learned symmetry discover the equation for the latent dynamics. We compare our approach with SINDY Autoencoder (Champion al. (219) without any symmetry and La SINDy, ie. learning the equations in the LaLiGAN latent space but without equivariance constraint.
Timestop
system and the pre- diction error in the high- high-dimensional input space using equations from different methods. The means and standard deviations (shaded area) of errors over random runs are plotted.
Figure 3:
Figure shows the prediction error on this system using the latent equations. We find that the discovery result of SINDy Autoencoder is unstable, and the simulation can quickly diverge. Meanwhile, both approaches with symmetry discover equations that accurately simulate the dynamics in the input space. Interestingly, even with similar latent spaces discovered by LaLiGAN, our approach that enforces the equivariance constraint can further reduce prediction error compared to non-constrained regression.
5.3 Symmetry Regularization
Table 2: Equation discovery statistics on the Lotka- Voltera system (15) at noise level o 99% and the glycolytic oscillator (16) at noise level σR 20%. The RMSE is scaled by x x10 for (15) and x10-2 for (16). The success probability is computed from runs for sparse regression-based algorithms and 10 runs for genetic mrogramming-based algorithms. The success probabilities of recovering individual equations (Eq. & Eq. 2) and simultaneously recovering both equations (All) are reported. The RMSE (all) refers the parameter estimation error over all runs. The RMSE (successful) refers to the parameter estimation error over successful runs, which is missing for algorithms with zero success probability.
Success prob.
(successful) Eg. 2
(all)
System
Method
Eq.
Eq. 2
All
Eq. 1
All
Eg.
Eg.
SINDy
0.40
0.64
0.24
1.01 (0.26)
0.56 (U.21)
(0.96)
(0.16)
4.01 (2.50)
(3.67)
WSINDY
0.22
0.32 (0.23)
(0.16)
(13.28)
(18.81)
EquivSINDy--
0.54
0.58
0.36
1.00 (0.21)
0.45 (0.20)
0.79 (0.15)
3.16 (2.46)
3.83 (4.01)
L-V
0.0
(0.89)
N/A
(0.89)
(3.20)
D-CODE
0.0
0.0
0.0
N/A
N/A
N/A
10:38 (0.
9.08 (1.88)
EquivGP-r
1.0
0.8
2.43 (1.39)
0.51 (0.98)
1.58 (1.62)
2.43 (1.39)
1.76 (0.31)
SINDY
0.30
0.36
0.14
0.87 (0.15)
0,32 (0.13)
0.67 (0.10)
15.86 (17.52)
12.71 (18.32)
Glycolytic
WSINDy
0.06
0.14
0.04
0.11 (0.10)
0.59 (0.23)
0.34 (0.04)
2.3e3 (2.9e3)
2.1e3 (2.7e3)
Oscillator
EquivSINDy--
0.40
0.70
0.28
0.92 (0.22)
0.30 (0.13)
0.71 (0.16)
9.97 (8.07)
7,29 (12.72)
16
0.0
0.0
0.0
NIA
NIA
N/A
7.57 (5.46)
14.85 (23.50)
D-CODE
0.0
0.4
N/A
0.27 (0.
N/A
14.37 (15.66)
8.90 (16.84)
EquivGP-r
0.1
0.6
0.1
0.67 (N/A)
0.22 (0.38)
0.40 (N(A)
8.59 (12.10)
3.43 (6.69)
We demonstrate the use of symmetry regularization in the following ODE systems. These systems do not possess any linear symmetries for solving the constraint explicitly.
Lotka- Lotka-Volterra System describes the interaction between predator and prey population. We consider its canonical form as in Yang et al. (2023a), where the state is given by the logarithm population densities of the prey and the predator.
2
22
3
(15)
+ 2L

---

Glycolytic Oscillator (Sel Kov 1968) is a biochemical system of two ODEs with complex interac- tions in cubic terms. We use the same constants in the equations as Qian et 2022).
#1 0.75 0.1.2 x102
(16)
22 0.1z1 22 + 2125
Since we assume no knowledge of the symmetries here, we first run LaLiGAN (Yang et al. 2023a discover symmetry from data for each system. We then extract the group action and the induced infinitesimal action from LaLiGAN and use them to compute symmetry regularization. Table 2 shows the results of different discovery algorithms with and without symmetry. Comparison is made within each class of methods, i.e. sparse regression and genetic programming. For genetic programming, we use an alternate form of symmetry regularization (37). which will be discussed in Appendix C.2
Informed by symmetries through regularization terms, equation discovery algorithms have a much higher success probability of discovering the correct equation forms. Also, models with symmetry regularization generally have the lowest parameter estimation emrors averaging over random experiments.
we provide additional results show that the equations discovered by our method also achieve lower prediction error. We also compare the different regularization options based on
In Appendix C.2
Theorem 3.3
Conclusion
We propose to incorporate symmetries into equation discovery algorithms. Depending on whether the action of symmetry is linear, we develop pipeline (Figure 1 with various techniques, c.g. solving the equivariance constraint and using different forms of symmetry regularization, to perform symmetry-inie -inforrr- equation discovery. Experimental results show that our method can reduce the search space of equations and recover true equations from noisy data.
Our methodology applies time-independent symmetries for autonomous ODEs. These symmetries have relatively simple forms and can be either obtained as prior knowledge or learned by existing algorithms. In theory, method be generalized to other settings, such time-dependent symmetries and non-autonomous ODEs. We provide a more general problem definition and discuss these potential generalizations Appendix [A.]] However, in these problems, less likely that we know the symmetries before discovering the equations. Extending our symmetry-intone methodology to these scenarios would require further breakthroughs in learning symmetries for general differential equation systems, which could be an interesting future research direction.
Acknowledgement
This work supported in part by the U.S. Army Research Office under Army-ECASE award W91 INF-22333-10023 the U.S. Department Of Energy, Office of Science, IARPA HAYSTAC Program, CDC-RFA-FT-23-00 DARPA AIE FoundSci, DARPA YFA, NSF Grants #2205093, #2146343 and #2134274.
References
Abou-Ismail, A. Compartmental models of the covid-19 pandemic for physicians and physician- scientists. SN comprehensive clinical medicine, 2(7):852-858, 2020.
Agapitos, A., Brabazon, A., and O' O'Neill, M. Controlling overfitting symbolic regression based on bias/variance error decomposition. In International Conference on Parallel Problem Solving from Nature, pp. 438-447, Springer, 2012.
Akhound-Sadegh, T., Perreault-Levasseue L., Brandstetter, Welling, M., and Ravanbakhsh, Lie point symmetry and physics informed networks. arXiv preprint arXiv:2311.00 1,44993, 2023.
Baddoo, P. J., Herrmann, B., McKeon, B. J., Kutz, J. N., and Brunton, S. I.. Physics-informed dynamic mode decomposition (pidmd), 2021.
10

---

Bakarji, J., Callaham, J., Brunton, S. L., and Kutz, J. N. Dimensionally consistent learning with buckingham pi. Nature Computational Science, 2:834 844, 12 2022. ISSN 2662-8457 doi: 00.1038864358800220000333
URI. https:///// http::://0oo..o/111111333333 1001:010000000/55000000000003330
Benton, G., Finzi, M., Izmailov, P., and Wilson, A. G. Learning invariances in neural networks from training data. Advances in neural information processing systems, 33: 7600-177666, 2020.
Brandstetter, J., Welling, M., and Worrall, D. E. Lie point symmetry data augmentation for neural pde solvers. In International Conference Machine Learning, pp- 2241-2250 PMLR, 2022.
Bronstein, M. M., Bruna, J., Cohen, T., and Veličković, P. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2104. 13478, 2021.
Brunton, S. L.,, Proctor, J. L.,, and Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(155:3932-3997, 2016. doi: 10, 10073//aa..111334411111 abs/10. 1011100.7//ppass... 1tw.::///p.....1.11111111 1517384113
URL BRLLhtts:://wer pnas org/doi/
Champion, K., Lusch, B., Kutz. J. N., and Brunton, S. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 166(4))22445-222 22451, 2019. doi: 10.1073/pnas.. 1903/99911000 URL https: pnas pnas. orgg//i/////ooae /0obs///.... 1073/ pnas 1906995116
Cohen, T. and Welling, M. Group equivariant convolutional networks. In International conference on machine learning, pp. 2990-2999, PMLR, 2016.
Cranmer, M. Interpretable machine learning for science with pysr and symbolicregression..J jl, 2023.
Cranmer, M., Sanchez Gonzalez, A., Battaglia, P., Xu, R., Cranmer, K. Spergel, D., and Ho, S. Discovering symbolic models from deep learning with inductive biases. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), Advances in Neural Information Processing Systems, volume 33, pp- Inc., 2020. URL htttps:// httpeeeeiiinnn.n.eni...re.eeeepp neurips cct.ccccppppee ffl_fiiesss wttps//ppppppp/iiiiiiiee 20FF2617777602226222226623338883 Paper Sraaer.pdi. pdf
17442. Curran Associates,
Cranmer, M. D., Xu, R., Battaglia, P., and Ho, S. Learning symbolic physics with graph networks, 2019.
Dehmamy, N., Walters, R., Liu, Y.. Wang. D., and Yu, R. Automatic symmetry discovery with lie algebra convolutional network. Advances in Neural Information Processing Systems, 34: 2021.
Finzi, M., Stanton, S., Izmailov, P., and Wilson, A. G. Generalizing convolutional neural networks for equivariance lie groups on arbitrary continuous data, 2020.
Finzi, M., Welling, M., and Wilson, A. G. A practical method for constructing equivariant multilayer perceptrons for arbitrary matrix groups. In International Conference Machine Learning, pp. 3318-332 3328. PMLR, 2021.
Gaucel, S., Keijzer, M., Lutton, E., and Tonda, A. Learning dynamical systems using standard symbolic regression. In Nicolau, M., Krawiec, K., Heywood, M. Castelli, M., Garcfa-Sánchez, P., Merelo, J. J., Rivas Santos, V. M., and Sim, K. (eds.), Genetic Programming, pp. 25-36, Berlin, Heidelberg, 2014. Springer Berlin Heidelberg.
Guan, Y., Brunton, S. L.. and Novosselov, I. Sparse nonlinear models of chaotic eletroroccctttooooveetttooo Soc. open sci., 8:202367, 2021. doi: 10. 20.10097/ URL http:://doi.oorg/0... r808. 202367
Huh, I., Yang. E., Hwang, S. J., and Shin, Time-reversal symmetric ode network. Advances in Neural Information Processing Systems, 19:16611-1 2020.
Kaba, S.-O., Mondal, A. K., Zhang, Y., Bengio, Y., and Ravanbakhsh, S. Equivariance with learned canonicalization functions. In International Conference on Machine Learning, pp. 15546-15566. PMLR, 2023.

---

Lee, K., Trask, N., and Stinis, P. Structure-prect preserving sparse identification of nonlinear dynamics for data-driven modeling. In Dong, B.. Li. Q.. Wang, L., and Xu, Z.-Q. J. (eds.). Proceedings of Mathematical and Scientific Machine Learning, volume 190 of Procedingss of Machine Learning mmlr press/
Research, pp. 65 -0. PMIR, 15-17 Aug 2022. URI. v190/1ee22a. html
Liu, Z. and Tegmark, M. Machine learning hidden symmetries. Physical Review Letters, 128(18): 180201, 2022.
Loiseau, J.-C. and Brunton, S. L. Constrained sparse galerkin regression. Journal of Fluid Mechanics, 838:42-67, 2018. doi: 10.1017/jfm.2117. 20fffm............
Messenger, D. A. and Bortz, D. M. Weak sindy: Galerkin-based data-driven model selection. Multiscale Modeling & Simulation, 19(3):1474-149) 2021.
Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., et al. Sympy: symbolic computing in python. PeerJ Computer Science, 3:3 2017.
Mialon, G., Garrido, Q.. Lawrence, H., Rehman, D., LeCun, Y., and Kiani, B. Self-supervised learning with lie symmetries for partial differential equations. Advances in Neural Information Processing Systems, (6:28973-29008 2023.
Nocedal, J. Updating quasi-newton matrices with limited storage. Mathematics of computation, 35 (151):773-782, 1980.
Olver, Applications of Lie groups to differential equations, volume 107. Springer Science & Business Media, 1993.
Otto, S. E., Zolman, N., Kutz, N., and Brunton, S. L. A unified framework to enforce, discover, and promote symmetry in machine learning, 2023.
Qian, Z., Kacprzyk, K., and van der Schaar, M. D-code: Discovering closed-form odes from observed trajectories. In International Conference on Learning Representations, 2022.
Ridderbusch, S., Offen, C., Ober-Blöbaum, S., and Goulart, Learning ode models with qualitative structure using gaussian processes. In 2021 60th IEEE Conference on Decision and Control (CDC), pp. 289222896999 2296, 2021. doi: 10. N0.1109/CDC45484..22 202.2.222..22.
Sel'Kov, Kov, H. Self-oscillations in glycolysis 1. simple kinetic model. European Journal of Biochem- istry, 4(1):79-86, 1968.
Sober, E. The principle of parsimony. The British Journal for the Philosophy of Science, 32(2): 145-156, 1981.
Udrescu, S.-M. and Tegmark, M. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16)eaay2631, 2020.
der Ouderaa, T., Immer, A., and van der Wilk, M. Learning layer-wise equivariances automatically using gradients. Advances in Neural Information Processing Systems, 36, 2024.
Wang, R., Walters, R., and Yu, R. Incorporating symmetry into deep dynamics models for improved generalization. In International Conference Learning Representations, 2021.
Weiler, M. and Cesa, G. General (2)-equivariant steerable cnns. Advances in Neural Information Processing Systems, 32, 2019.
Xie, X., Samaei, A., Guo, J., Liu, W. K., and Gan, Z. Data-driven discovery of dimensionless numbers and governing laws from scarce measurements. Nature Communications, 13()):762,, 2022. doi: 00.10884441467--2225555533
URL https:://oot org/ /0. 041467744167-0000 -22-35084--
Yang. J., Dehmamy, N., Walters, R., and Yu, R. Latent space symmetry discovery. arXiv preprint arxiv:2310.00105. 2023a.
Yang, J., Walters, R., Dehmamy, N., and Yu, R. Generative adversarial symmetry discovery. In International Conference Machine Learning, PMLR, 2023b.
12

---

A Preliminary on Lie Point Symmetry of Differential Equations
We use the theory of Lie groups to describe the symmetry transformations of ODE solutions. This section provides a brief overview some fundamental concepts of Lie group and its application to differential equations.
A.1 Lie Groups
Lie groups. A Lie group is smooth manifold with compatible group structure. We use Lie groups to describe continuous transformations. For example, all n x invertible matrices with real entries form the general linear group GL( IR) with matrix multiplication as the group operation.
The Lie algebra g == TeG defined as the tangent space at the identity element e of the group G. For example, The Lie algebra of general linear group GL(nn .) consists of all real-valued matrices of size n X n. We can use a vector space basis {v) E g} to describe any Lie algebra element as == Σ 6 where of the group. We can map the Lie algebra elements to group elements (in the identity component) via the exponential map, exp: g G. For matrix Lie groups, matrix exponential is such map.
R. Lie algebra can be thought of as the space of infinitesimal transformations
Group actions. A group action defines how abstract group elements transform other objects in a space V. The group action is a function I: G V
which maps the identity e to identity transformation, i.e. r(e, v) = and is compatible with group element composition, i.c. T(g.. (gg,, ")) = [gg1g2, v). We can also think of the group action as mapping group element to a function V V and write r == [gg,) -). The group action I induces an infinitesimal action y of its Lie algebra given by 7o (xx)) = d (I exp(ev) (z))) For a smooth manifold V, the infinitesimal action can be interpreted as vector field on that generates the flow coinciding with the group action of exp(ev). Given the coordinate (#1, .........a Fz), the vector field can be written Yo Σ, v'O,. where 0 Ss; and o² V R are the components of %o.
A group action is called representation if V is vector space and 9 is a linear function on V. When V = R", the representation r GL((n; R) is an invertible n x n matrix that transforms V E V by matrix multiplication. Correspondingly, the induced Lie algebra representation 7, E gl(n, R) is also an n x n matrix. the following discussion, we will refer the symmetries with linear actions as linear symmetries and those with nonlinear actions as nonlinear symmetries.
Notations. In this paper, we mostly consider the group action and the corresponding infinitesimal action on the phase space of the ODE, X. For simplicity, use g E G to refer to both the group element and the group action I' and write == I (æ). Similarly, we use E to refer to both the Lie algebra element and its infinitesimal action and write v(xr) = (x)).
A.2 Lie Point Symmetry of ODEs
A.2.1 ODE
In this work, we are mainly concerned with the autonomous ODEs in the form II. However, to introduce the general notion of Lie point symmetry of ODEs, it is helpful to consider a more generalized setup as follows:
dx = htt,x:)
(17)
dt
where é E T C IR denotes time, i.c. the only independent variable in ODEs, and x E X C Rd denotes the dependent variables.
To formalize the concept of symmetries, we take a geometric perspective on the differential equations (7). The independent and the dependent variables form a total space order jet space of the total space as J T x X x X1. Its coordinates, (1, (1,,1,..., In, 21 ..., rnn) E J. represent the independent and the dependent variables as well as the time derivatives of the dependent variables. This is also known as the first-order prolongation.
T x X. We define the first-
13

---

Then, we can represent the ODE T7 through the map A : J R" with components A; F/ (t, x). function ætt) solution of (17) its prolongation satisfies
af,, 0, Vt
(18)
A.2.2 Symmetry
symmetry of the system A transforms one of its solutions to another. More specifically, we need to define the group action on point in the total space E = 'I' x X and its induced action on function X.
We express the action of group element g in the total space as
g (f,z)= = (f(t, (t,,
(19)
where i and & are functions over T X.
To define the induced group action on the function f f T X, we first identify f with its graph
7f f()) t E Ωt) C T x X
(20)
where Qr C T is the time domain. The graph of the function can be transformed by as
g7f {ēē, x) g x)) (t. a²) E 7
(21)
The set gyy may not be the graph of another single- valued function. However, by choosing a suitable domain S2,, we ensure that the group elements near the identity transform the original graph into a graph of some other single-valued function æ = f(i) with the transformed graph V/ = -. In other words, the transformed function f g is defined by the transformation of the function graph.
Definition A.1 (Def 2.23, Olver 1993)) A symmetry group of A (88) is local group of transfor-
mations acting on an open subset of the total space I x X that whenever æ = f(t) is a solution of A, and whenever g f is defined for g G, then = (g (t) is also solution of the system.
A.3 The Infinitesimal Criterion
Let be vector field the total space T X with corresponding one-parameter subgroup exp(ev). The vector field describes the infinitesimal transformations to the independent and the dependent variables. It can be written as
vlt, x) [1, + Σ x1, 0 Fe,
(22)
The vector field also has its induced intinitesimal action on the jet space T x X x X2 described by the prolonged vector field:
(1) = = Σ
(2) (1,2,
(23)
where of be computed by the prolongation formula (Theorem 2.36, Olver (1993)
o (1)
=D,(0;
E.ril + Ex=!'
(24)
+ > x DΦ: DE dx, Ot +
ag dxx
x
Σ x
(25)
Ot
We denote (((, = 10 (1, ..., aa(1, x² the vector field components in X.
As a special case, for the time-independent symmetries we considered previously, { = and ott, æ) = ((æ), we have
0ø: dx, x/
= Σ
(26)
(1) Another special case is time translation. When & and Φ = 0, we have 0i
0.
The following infinitesimal criterion gives the symmetry condition in terms of the vector field (and its prolongation).
14

---

Theorem A.2 (The infinitesimal criterion. Thm 2.31, Olver (993)). Let A be a system of ODEs of
maximal rank defined over M CT X, X C RR. local group of transformations G acting on M is a symmetry group of the system A iff for every infinitesimal genrator v of G, (0) 0,
i= 1,.., d, whenever A 00.
Once we have the symmetry represented by the vector field v, we can derive constraints on the equation form from the above theorem.
Generally, when the vector field takes the general form in 22 Theorem A.2 is equivalent to
[(1, x)
(t) + J=(x)h(x)) (G. h()))) Jn (t, a))
(27)
((t,
For more comprehensive understanding of the application of Lie groups in differential equations. we refer the reader to Olver (993).
B Proofs
B.1 Lie Point Symmetry and Equivariance
Proposition B.1. Let G be group that acts on X, the phase space of the ODE ① G is symmetry group of the ODE ① in terms of Definition 1 if and only if for any T T, the flow map fr is equivariant to the G-action on X.
Proof. The proof follows from applying various definitions. If G is a symmetry group of (T. then for any solution = xt)), I = g is still solution to ( By definition of the flow map fr. we have (wo) Z.. where 3g = g (t; xo) and w, == g. z(t g 7), for any 2 E X. That is. f- (mo)) i.e. f- is G-equivariant.
fg To) = g +
On the other hand, for any solution the ODE ① xxt) where ==()) Fn, consider the transformed function क = x. Because fr G-equivariant, we have
#(t T) xtt + g fr (x((t) (ā(()), Vr
(28)
Taking the derivative w.r.t at t = 0, we have dt a(t) h()(t)), which means z is also solution to the ODE TT..
Theorem B.2. Let G be time-ine inneieddeeendent symmetry group of the ODE Let u E be an element in the Lie algebra g of the group Consider the flow map f defined (2) for a fixed time interval. For all x E XX and g = exp(ev) G. the following equations hold:
f(g m) 9 f(x) 0;
(29)
Jy(æ))u(x)) -f
v((())) - 0;
(30)
(g(x)((x) - h(g-2 x) = O;
(31)
Jn (n)))((z)) - (x)))((( = 0.
(32)
Proof. From Proposition B.1 the flow map is equivariant G. so (2) follows directly.
Consider a group element in the identity component that can be written as = exp(ev), v E g.
Taking the derivative with respect to on both sides of 29 at e = 0, we have
f(g x)]
gg f(z)
de
==0 de
==0
[f(g x)]
== v(f(z))
de
()((ep((0)(=)))(()) v v((()) = ()(())( = v(f())
and 32 can be proved similarly. Let (xo) to hxx dt (where xx(t) = 2)) be the flow map given the time interval r. Note that of (x) (x) hz). Since (2) is true for any т, we can take the derivative with respect
T 0 and obtain
Strictly speaking, the reverse only holds when A satisfy certain additional local solvability conditions.
15

---

Finally, we take the derivative with respect to on 30 at 0. We have
2 d
Jr. (ar)
r=
dry dr (=
dr
[Jn(z))
and
d
dr v(ff(())) J, -J((fro(()))hia = (ax)a)((
which give us 32
Remark B.3. 32) can also be obtained from the infinitesimal criterion (Theorem A.2 One can verify that when we assume time-independent symmetries and thus the prolonged vector fields given by (26). the symmetry condition in Theorem A.2 is equivalent
Σ doi dxy
Ōh
n'j E dxy
Dj.
i i= I=,...,d d
(33)
or, in its matrix form, (=)))((()) = JJ (x)) (z). where v(x) =
0a(x))
Remark B.4. (30) is similar to results in the existing literature, e.g. Section 4.2 in Finzi et al. (2021) and Theorem Otto et al. (2023). (30) is a generalization of these results which does not require the linearity of the function or the infinitesimal action
B.2 The Symbolic Map Me
In Section
4.1
we have defined the symbolic map Me (R"-R") + R*xp, which maps each
component of function with multivariate output to coordinate the function space spanned by O. This map is introduced to account for the fact that the constraint Proposition 4.1 holds for any We apply this map to search for equal functions by their equal coordinates in the function space.
As an example, let the input dimension be n = and θ the set of all polynomials up to second order, i.. (z), x2) [L, 1 32, z1, r122,
xz[T" 1222, Define 22)) = [xž] x²)" and fa(21, x2) z/.
0
0
0 1
We have Me (f11
0 and Me(fz) undefined because the cubic term cannot
be expressed in terms of linear combination of functions in O.
Thus, for a given infinitesimal action v, we need to ensure that an appropriate function library θ is chosen for SINDy that it can express the equivariant functions to the given action. For linear action v(x) L,æ, we prove that the polynomial functions satisfy the requirement.
Proposition B.5. Let (xx) be the set of all polynomial functions up to degree q (g E EZ). Then, the Lo Me(Je())Lo()) is well-defined.
components of Je) L E Rp can be written as linear combinations of the terms in O(x), so that
Proof. Consider the function library (xx) that includes all polynomial terms up to degree p 2 1. The Jacobian Je (x: captures the partial derivatives of each component of (zx), which are polynomials up to degree polynomial in terms of æ. Hence, the product Je(x) L. yields function vector whose components are sums of products of two sets of polynomials: those of degree up to q - from the Jacobian, and from the linear transformation. The highest possible degree of any term in the degree one.
1. The linear transformation L, maps æ to L,æ, and therefore, is a first-degree
those of degree
product is q, as is the sum polynomial of degree q and a polynomial
Since (x)) includes all polynomials up to degree p. it includes all the terms generated by the product Jel (a:) T. Therefore, the components of Jol (z) T, E TRP can be written as linear combinations of the terms in (xx)
16

---

C Supplementary Experimental Results
C.1 Equivariant SINDy for Linear Symmetries
C.I.1 Supplementary Results for Section 5.1
In addition to Table in the main text, Table B the parameter estimation error and the long-term prediction error based on the discovered equations. For the parameter estimation error, we report the RMSE over successful runs and all runs. For the long-term prediction error, we select three temporal checkpoints and report the average prediction error of the equations at these checkpoints. The predictions from our discovered equations are much more accurate than equations from other methods. While our EquivSINDy-c achieves the best success probability and parameter estimation eπrоr over all runs, it is interesting to observe that Weak SINDy (WSINDy) has very accurate estimations of the constant parameters once it discovers the correct equation form (occasionally). We conjecture that this is because using the weak form of ODE averages out the white noise over long time interval. Its low success probability suggests that the objective based on the weak form may introduce some difficulties in eliminating the incorrect terms. However, once it is given the correct form of the equation, parameter estimation can become much more accurate compared to other methods that use L2 error at each individual data point with high measurement noise.
Table 3: Equation discovery statistics on the damped oscillator 13 at noise level dR 20% and the growth system (14) OR = 5%. The success probability is computed from 100 runs for each algorithm. The success probabilities of recovering individual equations (Eq. & Eq.2) and simultaneously recovering both equations (All) are reported. The RMSE (all) refers to the parameter estimation error over all runs. The RMSE (successful) refers to the parameter estimation error over successful runs, which is missing for algorithms with zero success probability.
Success proh.
RMSE (successful) (x 10
1
RMSE (all) (x 10 1
System
Method
Eq.
Eq.
All
Eq
Ey. 2
All
Eq. 1
Eq 2
GP
0.00
0.00
0.23 (0.II)
NIA
0.71 (..00
0,37
D-CODE
0.00
0.00
0.00
N/A
NIA
N/A
7.59 (9.35)
7.53 (5.44)
Oscillator
SINDY
0.35
0.15
(0.06)
(0.09)
(0.04)
0.50 (0:22)
0.48 (0.21)
WSINDY
0.06
0.07
0.00
0.10 (0.05)
0.09 (0.05)
N/A
1.17 (0.56)
1.07 (0.57)
EquivSINDy-c
0.93
0.97
0.90
0.37 (0.07)
0,37 (0.08)
0.37 (0.08)
0.37 (0.08)
0.37 (0.08)
GP
0.00
0.00
(0.14)
2.13 (0.00)
0.16
D DODE
0.00
0.65
0.00
NIA
1.16 (0.28)
N/A
2.17 (0.26)
0.87 (0.55)
Growth
SINDY
0.26
0.03
(0.05)
(0.02)
(0.51)
0.50 (0.46)
WSINDY
0.00
0.00
0.00
N/A
N/A
N/A
2.45 (0.14)
4.77 (3.34)
EquivSINDy--
1.00
1.00
1.00
0.19 (0.10)
0.08 (0.06)
0.15 (0.07)
0.19 (0.10)
0.08 (0.06)
Long-term orediction error for damped oscillator
Long -term prediction error for growth system
0,00
Eqqivs INJyCc
EquivSINDy-c
SINDY
0.025
SINDy
GF
0,020
0,015
0.3
MSE
02-
0.000
-0.003
Figure 4: Long-term prediction error (MSE) against simulation time. The error curves are averaged over all discovered equations and random initial conditions in the test dataset. The shaded area indicates the standard deviation of prediction error at each timestep. Our EquivSINDy-e has the slowest error growth. Some algorithms, c.g. genetic programming and Weak SINDy, are not included in the plot when the simulation quickly diverges and the error grows to infinity.
Figure plots the error growth curve when using the discovered ODEs to simulate future dynamics. Some algorithms, e.g. genetic programming and Weak SINDy, are not included in the plot when
17

---

the simulation quickly diverges and the error grows to infinity. We observe that the prediction error grows slowest when using the discovered equations from our method, EquivSINDy-c.
C.1.2 A Higher-Dimmensonall System
Our method can also be applied to higher-dimensional dynamical systems with symmetries. To demonstrate this, we consider the SEIR equations Abou-Issaaa Ismaill (2020) which model the evolution of pandemics by the number susceptible (S), exposed (E), infected (D, and recovered (R) individuals. It is dimenssional ODE system with quadratic terms. We the following equations the target of equation discovery:
$ = 0.15 0.651
É - 0.6S1
(34)
i= E - -0.5I
R = 00.55 + 0.5L
Compared to the original formulation, we add two constant terms to the first and the fourth equation to account for the re-infection of patients due to waning immunity, mutation of the pathogen, etc. During data generation, we add moderate 5% noise all trajectories.
We assume an intuitive symmetry of scaling the number of recovered individuals proportional to the total population: -S E I + R)0g. We solve the constraint with regard this symmetry as Figure shows the reduced parameter space (from 60D to 34D).
in Section
Figure 5: In the SEIR model, we search the equations of variables with up to quadratic terms, leading 4x 15=60D parameter space. The symmetry (S I R)On reduces it to 34D.
4 the discovery results of SINDy and our EquivSINDy-C terms of success probabilities and parameter estimation eтrоr when successful. It turns out that SINDy struggles to find the correct equations and never manages to correctly identify all the equations simultaneously. In comparison, by reducing the parameter space through symmetry constraint, our method greatly increases the chance of successful discovery.
Table 4: Equation discovery SEIR epidemic model (4D). Results evaluated 50 random runs.
Success prob.
RMSE (successful)
Dynamics
Method
Individual eqs
Joint
Individual eqs
Joint
SINDy 0.02/0.02/0.20//0 10
0.00
1.57/0.3//.42/11..
SEIR
EquivSINDy-c
0.26/0.98/0.84/066
0.14
1.22/0.55/1 8111166666666666 1.6
1.57 (0.47)
C.2 Symmetry Regularization
In addition to Table 2 the main text, Figure 6 and Ssow the error growth when using the discovered ODEs to simulate future dynamics. Some algorithms are not included in the plot because the simulation from the learned equation quickly diverges and the error grows infinity. For the Lotka- Volte-Vo system in Figure 6 we observe that the prediction error grows slowest when using the
18

---

Long-term prediction error for Lott---attrr system
EquivSiNDy
SINDy
34+
EquixCil>r
48
0.00
0.35
3.50
6.75
3.00
1.25
3.50
3.75
2.00
Time
Figure 6: Long-term prediction error (MSE) on Lotka- Voltera system. The error curves are averaged over all discovered equations and random initial conditions in the test dataset. The shaded area indicates the standard deviation of prediction error at each timestep. Our methods based on symmetry regularization have slower error growth (visualized in solid lines).
discovered equations from methods based symmetry regularizations. Among all the methods, genetic programming with symmetry regularization has the most accurate prediction.
Figure 7 shows the results on the Sel' Kov glycolytic oscillator. Our symmetry regularization improves prediction accuracy over both of the base methods (SINDy and GP). However, D-CODE can achieve the best prediction accuracy in long horizon, even though its discovered equations may not have the correct form.
Besides, it is noteworthy that genetic programming with symmetry regularization leads to fluctuating error curve (after averaging). A potential reason that may discover some equations that are far from accurate but still lead to periodic trajectory, similar to the true trajectory of the Sel' Kov glycolytic oscillator. The trajectory may have large "eccentricity", so the error grows fast during a certain time and decreases rapidly afterward.
Prediction for glycolytic oscillator
Eqyi
SINDy
6.8
EquivGPrr
D-CODE
6.4
6.2
10.0
17.5
Time
Figure Long-term prediction error (MSE) on glycolytic oscillator. Left: results for sparse pegression-based methods. Right: genetic programming based methods. Our symmetry regularization improves prediction accuracy over the base methods. D-CODE can achieve the best prediction accuracy in long horizon, even though its discovered equations may not have the correct form.
Besides, we also visualize the simulated trajectories some learned equations in Figure S This provides more intuitive explanation of the error growth. For example, the exploding variance in the long-term prediction error in the Lotka- Volterra system can be seen in Figure Sa where predicted trajectories from most of the learned equations stay close to ground truth, but few equations lead to diverging predictions due to the incorrect modeling of the exponential terms.
Regularization Options. mentioned Section we can introduce different forms of symme- try regularizations based on the four equations in Theorem 3.3 listed as follows. The subscripts
19

---

Simulated trajectories from EquivSINDy.r
trajectories EquivSINDy
(a) Lotka- Volterra system
(b) Glycolytic oscillator
Figure 8: Simulated trajectories of the learned equations from EquivSINDy-r. In each figure, we sample learned equations and start the simulation with the same initial value. The solid lines of different colors correspond to different discovered equations.
short for Infinitesimal/Fiiti Group actions and Infinitssimll/Finiettiie flow map of the ODE.
|| -
LIGFE
Ex
(35)
|| (x))x()))
vEB(g)
(g -) -
ff f2)
PROFE = Ex
Σ
(36)
| (g x) - (2)
vE
[1((2)((z) h1g x x)1
LFGIE = Ex
Σ
(37)
())()) 2
VE (()
[J((())(()) - J J=(2) [J(())((>)|2
LIGIE = Ex
(38)
((
where y exp(er) is computed for some fixed scale parameter e. In other words, the summation is always performed over the learned infinitesimal generators (u E (g)) and the group elements g in the above formulas are dependent v.
We also note that the denominator of (36) has a somewhat different form from others, i.c. the difference between fgæ) and fxx). If we f(gx) or fxx instead, the scale of the denominator would also depend on æ. As result, regions closer to the origin in the state space have larger weights (smaller denominators) during loss computation. It's helpful to consider what would happen in the limit of + 0,7 + 0. The formulas inside the norms in the numerator and denominator in (35) both become (cTr), while (gx) f (exp(cv) does not match the numerator.
+ f(x + O())) = + O(c) O((r), which
While the corresponding constraints in Theorem 3.3 of these losses are all necessary conditions for symmetry, these loss terms have different implementations that may bring certain advantages or visadvantages depending the application, e.g. the representation of symmetry and the algorithm for equation discovery.
for training. Empirically, we find that (35 36 37) improve the discovery performance similarly while (38) leads to a performance drop. Detailed results are available in Table[5] We conjecture this is because we train a neural network to model the finite symmetry transformations (discussed in Section 4.3 Computing the Jacobian J, therefore involves computing second-order derivatives of the network and amplifies the numerical error. On the other hand, we have exact knowledge about the symmetry of the system in terms of its infinitesimal generators
In practice, we use only one of the losses 35
36
38 becomes the most efficient way incorporate symmetry loss,
20

---

without the need of computing group elements based on the Lie algebra generators or integrating the learned ODE its flow map.
Table 5: Equation discovery statistics of sparse regression with different regularization options
35 36 37 38 Lotka- Volterra system 15 at noise level 0g 99% and Sel'Kov oscillator [6
at noise level 0R 20%.
System
Regularization
Success proh.
RMSE (successful)
RMSE (all)
Eq.
Eq.
All
Eq.
Fq.
Fq.
Eq. 2
IGFE
0.54
0.58
0.36
1.00 (0.21)
0.45 (0.20)
0.79 (0.15)
(2.46)
(..3(1
L-V 15
FGFE
0.54
01.6
0.34
0.99 (0.21)
0.45 (0.20)
0.79 (0.15)
3.29 (2.65)
4,00 (4.05)
FGIE
0.58
0.58
0.38
0.99 (0.20)
0.44 (0.20)
0.79
2.90 (2.28)
3.80 (4.00)
IGIE
0.48
0.16
0.08
0.58 (0.31)
0.66 (0.17)
0,88 (0.06)
3,31 (2.67)
7.61 (3:42)
IGFE
0.40
0.70
0.28
0.92
0.30 (0.13)
0.71
9,97
7.29 (12.72)
Sel'Kov
PGFE
0.48
0.68
0.28
0.87 (0.22)
0.26 (0. 12)
0.63 (0.13)
8.38 (8.55)
7:97 (14.60)
TG
FGIE
0.48
0.54
0.28
0.99 19)
0.72 (0.15)
(0.08)
8.79 (8.94)
10.85 (17.50)
IGIE
0,50
0.48
0.22
0.94 (0.17)
0.89 (0.16)
0.91 (0.09)
8.25 (8.09)
14.56 (20.21)
Besides, B7p and (88) have specific advantage that any computation related to group actions is independent of the function h which is constantly updated throughout training. Thus, we can compute the (infinitesimal) group actions and their Jacobians in advance over the entire dataset. This provides better compatibility between the use of symmetry and existing tools for equation discovery. For example, in PySR framework (Cranmer) 2023), is difficult to compute group transformations g x based on an external neural network module during the equation search. We overcome this challenge by pre-computing the transformations and using the transformed data as auxiliary inputs.
Table 6 summarizes the comparison between these losses. While their corresponding constraints in Theorem B.3 are all necessary conditions for symmetry, they have different implementations that may bring certain advantages or disadvantager depending on the application, e.g. the representation of symmetry and the algorithm for equation discovery. For the experiments in Section 5.3 for sparse regression and LIGIE for genetic programming.
use LIGFE
Table 6: Comparison between different symmetry regularization losses whether they need to compute (finite) group elements, compute higher-order derivatives (for infinitesimal action), integrate learned equations, and whether they can pre-compute the symmetry transformations over the data.
Loss
335)
(36)
(37)
(38)
Compute group elements
No
Yes
Yes
No
Higher-order derivatives
No
No
No
Yes
Integrate learned equations
Yes
Yes
No
No
Pre-compute symmetry
No
No
Yes
Yes
C.3 Equivariant SINDy in Latent Space
Table 7: Samples of discovered equations from SINDyAE Champion et al. 20199. LaLiGAN + SINDY Yang et al. (023a) and our Equivariant SINDyAE.
Method
Equation Samples
21
+ 16
2: 0.4223
à 1.085, + 0.325
SINDYAF
& = 0.16s2 0.80s
& 0.88s, 0.37s + +085
2g 0.28z)
$1 0.10s1
0.97%
&1 0.9952
31 0.11s1
+ 0.9292
LALIGAN+SINDy
82 = -0.88s,
22 -0.8681
82 = -0.925, 0.1122
21 -0.8882
21 -0.9922
0.78a2
EquivSINDyAE
#2 = 0.9221
$2 = 0.83s1
22 = -1.04s1
In Table show some latent equations for the reaction-difflsion system in Section 5.2 discovered by different methods. Note that there is no single correct answer for the latent equations, because there might exist different latent spaces where the high-dimensional dynamics can be described accurately via different equations. Thus, the analysis in this subsection is mainly qualitative.
The discovered equations from SINDY Autoencoder (SINDyAE) contain numerous quadratic terms. While these equations may achieve a relative low SINDy error evaluated on individual timesteps, they fail to account for the periodic nature of the reaction-diffusion system. As result, if simulate
21

---

the system towards long horizon with these equations, the enror will accumulate rapidly, as shown in Figure E] Also, the discovered equations from multiple runs are essentially different, making difficult to select the best one to describe the dynamics.
Compared to SIN symmetry- hased methods such as LaLiGAN + SINDY and our equivariant SINDy Autoencoder can stably discover simple linear equations with low prediction error. Moreover, we observe that our method yields more consistent results than the non-equivariant approach (LaLi- GAN + SINDy). In LaLiGAN + SINDy, although the latent space has a linear symmetry action, the discovered equations are not constrained satisfy the symmetry. The resulting equations can take slightly different forms, as shown in T 7
On the other hand, our method constantly discovers equations in the form & -222 ž2 bā1, where the numerical relationship between and b depends on the discovered symmetries.
C.4 Comparison with Neural-ODE. Dreral-OD--Baeee Methods
We compare our method with the Reversll Symmetry ODE Network (TRS-ODEN) (Huh et al. 2020), which uses the discrete symmetry of time-reversal to improve dynamics learning. TRS-ODEN aims to learn the physical dynamics with black-box neural networks and adds a loss term encourage time-reversal symmetry. We evaluate TRS- ERS-ODEN as well as other baselines in Huuh et al. (2020). i.e. Neural ODE (ODEN) and Hamiltonian Neural ODE (HODEN), on the damped oscillator system (13) and the Lotka- Voltera system (15) in terms of the long-term prediction error. Other metrics, including success probability and parameter estimation error, are not compared because Neural ODE does not learn a symbolic equation.
Long term prediction error for damped ascillator
Long-tc m prediction error for Lottr-Voltera system
6.4
LquivSINDy
EquivSINDy
0.35
SINDY
SINDY
NODE
EquivGPrr
HNODE
TRS-NODE
6.20
NODE
0.2
0:15
6.20
0.30
5.0
7.5
10.0
12.5
15.0
175
0.3
0.70
1.50
3.25
3.50
Time
Time
Figure 9: Long-term prediction errors of our method & Neural-ODE-based methods on the irreversible damped oscillator (left) and the reversible Lotka- Volterra system (right).
Figure plots the long-term prediction error against the prediction horizon. In the damped oscillator, because it is irreversible and not conservative, TRS-ODEN and HODEN do not help. In the reversible Hamiltonian system of Lotka-Volterra equations, TRS-ODEN and HODEN improve upon the baseline Neural ODE, but have slightly higher errors than discovered equations with Lie symmetry regularization. Notably, in the Lotka- Voltera system, our method does not rely on any prior knowledge (e.g. time-reversal symmetry conservation of energy) but still achieves lower prediction error. To conclude, these results demonstrate that our method can be universally applied to various dynamical systems and achieve superior performance.
D Experiment Details
D.1 Data Generation
Damped Oscillator. We generate 50 trajectories for training, 10 for validation, and 10 for training.
with T = 100 timesteps sampled at a fixed rate A/ 0.2. The initial conditions are sampled by sampling the distance to origin r ~ 2((0.5, 2] and the polar angle 0 ~ [0, 2nr] and then converting to (x1, .22) = (r cos 0, r sin 0). The numerical integration performed with step size of 0.002. We apply white noise with noise level og 20% to each dimension of the observation.
22

---

Growth. We generate 100 trajectories for training. 20 for validation, and 20 for testing, with T = 100 timesteps sampled a fixed rate At 0.02. The initial conditions are sampled uniformly from [0.2, 1] x [0.2, I]. The numerical integration is performed with a step size of 0.002. Instead of applying white noise, we apply a multiplicative noise to this system, so the noise grows together with the scale of the state variables as the system evolves. Specifically, the observation is obtained by ã = 21 + E;), where 2; is the true state from numerical integration and e is sampled from Gaussian distribution with zero mean and standard deviation 0.05. We still refer to this as 5% noise level in the main paper. While the assumption of Gaussian process does not hold here, we find that Gaussian process smoothing can still improve the likelihood for the equation discovery algorithms to recover the correct equation. Therefore, just as in any other ODE systems considered in this work, we apply Gaussian process smoothing before running equation discovery methods.
Lotka-Volterra System. We generate 200 trajectories for training, 20 for validation, and 20 for testing, with T = 10000 timesteps sampled at fixed rate At 0.002. For initial conditions, we first sample the population densities (pI. P2) uniformly from [0, 1] x [o, 1]. Then they are converted to the logarithm population densities in the canonical form of the equation, i.e.
log(p1) and Ø2 log(p2). Also, we ensure that the Hamiltonian of the system H = exp(x1) - I1 + 1.333 exp(22) 0.667 under the sampled initial conditions falls within the range of [3, 4.5). We discard the initial conditions that do not meet this requirement. Then, the numerical integration is performed with step size of 0.002. We apply white noise with noise level oR 99% to each dimension of the observation.
Glycolytic Oscillator. We generate 10 trajectories for training, 2 for validation, and 2 for testing, with T 10000 timesteps sampled at a fixed rate A/- 0.002. The initial conditions are sampled uniformly from [0.5, 1] x [0.5, 1]. The numerical integration is performed with a step size of 0.002. We apply white noise with noise level oR = 20% to each dimension of the observation. We note that our experimental results (2022) on this system with the same noise level. This is because our initial conditions are sampled from a larger range which leads to more diverse trajectories. As result, the metrics such success probabilitics and parameter estimation errors appear worse than in Qian et al. (2022). However, we still ensure fair comparison we use the same generated dataset with configurations for all the methods involved in our experiment.
this system are somewhat different from the results in Qian et al.
Reaction-Diffusion System.
The system is governed by the PDE
=(1 (u²² + +z)) (u² + + 0.1(uzz + Myy).
- (u² + v²) )u + (1 - (u² + )²)) 0.1 + Uyg))
We use the code from SINDY Autoencoder + Champion et al. 2019) to generate the data, The
2D space of u and v is discretized into a 100 x 100 grid. The system is simulated from one initial value for T = 6000 timesteps with step size At 0.05. Then, we add random Gaussian noises with standard deviation 10 to each pixel and at each timestep. We use the timesteps t E [0, 4800) for training. For forecasting, we use the timestep t = 4800 as the initial value and simulate up to 600 timesteps with each method. The simulations are compared with the ground truth during t e (4800, 5400) to calculate the error.
D.2 Equation Discovery
D.2.1 Sparse Regression
For all algorithms based on sparse regression (SINDy), we use a function library ((x) that contains up to second-order polynomials, with an exception for Lotka- VolteraVolteraaaaa system where we also include the exponential terms. We apply sequential thresholding (Brunton et al. (016) to enforce parsimony in the equations. The threshold is set 0.05 for the damped oscillator and the growth system, 0.075 for the glycolytic oscillator, and 0.15 15 for the all methods including SINDy, WSINDy and EquivSINDy.
system. The same threshold is applied to
The original SINDy uses a least square objective and can be solved explicitly. However, when we introduce symmetry regularization into the training objective, the original solver is longer
23

---

applicable. For our equivariant models based on symmetry regularization, we use the L-BFGS algorithm (Nocedal 1980) for optimizing the SINDy parameters. We find L-BFGS more effective than stochastic gradient descent for these low-dimensional problems.
D.2.2 Genetic Programming
We incorporate our symmetry regularization 37 to an existing genetic programming package for discovering equations, PySR (Cranmer) 2023 The "GP" baseline in the tables also refers to the implementation of genetic programming with PySR. On the other hand, we run the experiments This codebase is developed based another genetic programming package, gplearn. We should note that the performance difference between our EquivGP and D-CODE can be partially caused by their different backends.
on D-CODE with their official codebase |
24

---

NeurIPS Paper Checklist
The checklist is designed to encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No], or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justifiicton right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation.
While "[Yes] " is generally preferable to "[No]", it is perfectly acceptable to answer "[No] provided proper justification is given (e.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used"). In general, answering "[No] or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We have listed our four main contributions in the introduction, (roughly)
corresponding to Section 3] Section 4.1 Section 4.2 to Section 4.3 and 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 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 as 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: Our paper is limited to the scope of autonomous ODEs. See Section 6 for discussion on limitations and potential generalizations.
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.
25

---

The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with 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 low or images are taken in low lighting. Or speech-to-tex: 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 faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an 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 of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Appendix A discusses the scope of this paper and the related assumptions. Proofs to the theoretic results are available in Appendix B
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 or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We provide details about experiments in Appendix D codebase available in the supplementary material.
We also make our
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a 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 provided or not.
26

---

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 a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it 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, other means that 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 to reproduce that algorithm. (b) If the contribution is primarily new 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 to 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.. to registered users), but it should be possible for other researchers have some path 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 faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We make our codebase, including data generation and experiments, available in the supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
for more details.
The authors should provide instructions 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. 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 in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
27

---

Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: Section 5is the experiment settings in general and Appendix D provides more details.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to 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
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Each experiment is repeated for 10~50 times. The success probability metric is an aggregation over the replications, so no error bar is available. For other metrics, report the mean and standard deviation over all the replicated runs.
Justification:
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, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
run 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).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but 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 to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures 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 execution) needed to reproduce the experiments?
Question:
Answer: [No]
We did not keep track of the execution time as most experiments are solving Justification: low-dimensional problems and can be done within minutes on single CPU or GPU.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
28

---

The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make 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
Answer: [Yes]
Justification: The authors have read the NeurIPS Code of Ethics and confirmed that this paper does not violate
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there 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: [NA]
Justification: This paper aims to discover governing equations in physical science and has minimal societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has socictal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly. harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional 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 to monitor how system lears from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: The data used in the paper are generated by numerical methods and have no safety risks.
29

---

Guidelines:
The answer NA means that the paper poses no 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 to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but encourage authors to take this into account and make best faith effort.
12. Licenses 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 paper has cited the original papers that produced the dataset and provided the URLs.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g.. website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, phcorswithcode thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors are encouraged to reach out 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: [Yes]
Justification:
We provide the documentation along with the codebase in the supplementary
material.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model 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 used.
At submission time, remember anonymize your assets (if applicable). You 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)?
30

---

Answer:
[NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing 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
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?
Question:
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor 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.
31


Reviews:
[{"review_id": "QTXdNmF027", "reviewer": "Reviewer_jBog", "summary": "This paper focuses on incorporating symmetries to equation discovery pipelines for ODEs. \n\nIf the governing equation has a known linear symmetry (its solutions are invariant with respect to a known linear action by a known Lie group), this paper derives a set of conditions the ODE needs to satisfy, and incorporates those conditions into a symbolic regression framework such as SINDy. \n\nWhen the symmetry is still known but not linear, the paper derives a form of regularization to promote that symmetry into the learning process. This procedure wasn't completely clear to me.\n\nWhen the symmetry is not known, the paper uses LaLiGAN, a symmetry discovery framework, to learn the corresponding regularization.", "questions": "What are the limitations of the approach in the known non-linear symmetries case?\n\nWhat are the limitation of the approach in the unknown symmetries case?", "limitations": "The authors could comment on the technical limitations of their work (e.g. dimensionality, ODE order, sample complexity, etc).", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is a well-written paper in an interesting topic: imposing symmetries in governing equation discovery. \n\nIt considers approaches from the simpler case (known linear symmetries) to the more general case (unknown non-linear symmetries), and it provides numerical evaluations for them.", "weaknesses": "The methodology behind the regularization approach and the LaLiGAN approach were not clearly explained in the paper. It would be good if the authors can provide a more detailed explanation.\n \nIn particular the technical approach for unknown non-linear symmetries (based on LaLiGAN) was not clear to me. I went to the recent LaLiGAN paper to see how it works and I found that the original LaLiGAN paper already proposes the approach for governing equation discovery. I suggest the authors explain how the approach they propose in this paper differs from the approach proposed in LaLiGAN. Moreover, the LaLiGAN paper does not have a github repository with their code, so the authors should be able to make sure that all their code and their dependencies can be made publicly available if the paper gets accepted.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback. We address key comments below.\n\n> On the connection and difference from LaLiGAN\n\n**Our approach has a completely different goal from LaLiGAN.** Our method aims to discover equations using symmetry as an inductive bias. LaLiGAN aims to discover unknown symmetry.\n\nConnection: Our method requires the knowledge of symmetry. We consider two scenarios: one where symmetry is given to us as a prior, the other where we do not know the symmetry. Only in the latter case, **we use LaLiGAN as a tool to discover the symmetry of dynamical systems.** Our method can also work with any other symmetry discovery techniques.\n\nIn terms of equation discovery, we propose equivariant models for this task given the symmetry, while LaLiGAN did not discuss equivariant models or any other new method for equation discovery.\n\nTo explain the differences more specifically:\n* In Sec 5.1, we solve the equivariance constraints for linear symmetries, such as rotation and scaling. This is one of our new contributions and has not been discussed in LaLiGAN.\n* In Sec 5.2, we use LaLiGAN to discover a latent space (symmetry discovery) and then discover an equivariant equation in the latent space. In the second step (equation discovery), we enforce the equivariance constraint on the equation, whereas LaLiGAN uses a non-equivariant model to discover the equation. In other words, if we view the latent dynamics as the data for equation discovery, LaLiGAN has **equivariant data + non-equivariant model**, while we have **equivariant data + equivariant model**. In Figure 3, we show that using such an equivariant model that matches the symmetry of the data could further improve accuracy.\n* In Sec 5.3, we use LaLiGAN to discover nonlinear symmetry (same as LaLiGAN). Then, we encourage our model to be (approximately) equivariant to this symmetry through regularization (Sec 4.2). LaLiGAN did not use an equivariant model to discover equations. Same as in the above bullet point, they used equivariant data (in the latent space) + non-equivariant model. Moreover, they only learned equations in terms of latent variables. Our method discovers equations in the original state space, which is more interpretable.\n\nWe will clarify these differences in the updated version of the paper.\n\n> Code availability\n\nIn the Supplementary Material, we have included our implementation of LaLiGAN. Our code is fully runnable without external non-public dependencies.\n\n> Limitations in the case of known nonlinear symmetry / unknown symmetry\n\nWhen a nonlinear symmetry is known, we can directly apply the symmetry regularization approach, as shown in Sec 4.2. This would be a relatively easy scenario compared to our experiments in Sec 5.3, because we don’t need to discover the symmetry. There might be limitations to the regularization approach itself. For example, choosing the regularization coefficient may be cumbersome; the optimization problem may become more difficult due to non-convexity.\n\nWhen the symmetry is unknown, we use symmetry discovery methods such as LaLiGAN to learn the symmetry. If the discovered symmetry is inaccurate, it would cause an additional source of error in computing our symmetry regularization. Thus, our method indeed relies on the accuracy of the learned symmetry.\n\n> Other technical limitations\n\nWe thank the reviewer for pointing out these potential limitations. We provide some brief comments below and will discuss these points in more detail in the revised paper.\n\nRegarding dimensionality, our current experiments mainly focus on low-dimensional problems. Symmetries in higher-dimensional systems can be more challenging to identify. However, once we know the symmetry, our method can be easily extended to higher-dimensional problems. Also, higher-order systems can be reduced to lower order by introducing the higher-order derivatives as new state variables.", "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": "5nRUPArW5e", "reviewer": "Reviewer_wQw7", "summary": "The paper proposes to leverage symmetry to guide the equation discovery process, compress the equation search space, and improve the accuracy and simplicity of the learned equations. Depending on the types of symmetries, the paper develops a pipeline for incorporating symmetry constraints into various equation discovery algorithms, including sparse regression and genetic programming", "questions": "1. Why set 5% and 20% noise respectively in the experiment in Section 5.1? What happens when the error scale is smaller?\n\n2. For equation discovery with unknown symmetry, if we assume that the symmetry is unknown in the first two examples, what will be the result? \n\n3. The paper only considers two-dimensional situations. What are the results for high-dimensional examples?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper clearly establishes a pipeline to use Lie point symmetries of ODEs to constrain the equation learning problem.\n\n2. The paper theoretically derives the criterion for symmetry of time-independent ODEs in terms of equivariance of the associated flow map. \n\n3. From the above mentioned criterion, the paper solves the constraint explicitly to compress the equation search space in sparse regression and promote parsimony in the learned equations.\n\n4. In experiments across many dynamical systems with substantial noise, the symmetry-informed approach achieves higher success rates in recovering the governing equations.", "weaknesses": "The Time-Reversal Symmetric Ordinary Differential Equation Network has been proposed [1]. However, this work does not cite the relevant article [1]. While the presentation may differ, the symmetries investigated in both studies are fundamentally the same. Although this work is based on SINDy and  the relevant article [1] is based on Neural ODE, this may not be an essential difference. Therefore, the reviewer recommends a detailed comparison of the two approaches, highlighting their similarities and differences. Additionally, incorporating comparisons into the experimental section would be beneficial.\n\n[1] Huh, I., Yang, E., Hwang, S. J., & Shin, J. (2020). Time-reversal symmetric ode network. Advances in Neural Information Processing Systems, 33, 19016-19027.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback and the recognition of our paper’s presentation and evaluation. We address the key comments below.\n\n> Related work: Time-reversal symmetric ODE Network (TRS-ODEN)\n\nThank you for bringing this to our attention. This work aims to learn the physical dynamics with black-box neural networks whereas our goal is to discover the *symbolic* form of the governing equation. Furthermore, the range of symmetries considered in these two works differs. We consider the Lie symmetries of the ODE systems, which cover the broad range of continuous symmetry transformations.TRS-ODEN focuses on the specific symmetry of time reversal, which is a discrete symmetry. Because of the distinct natures of these symmetries, the ways of incorporating them also differ.\n\nFor long-term prediction tasks, TRS-ODEN can be used as a baseline. In **Figure 5** (supplementary pdf), we tested TRS-ODEN and other NeuralODE-based methods on two tasks in our paper. In the damped oscillator, because it is irreversible and not conservative, TRS-ODEN and HODEN do not help. In the reversible Hamiltonian system of Lotka-Volterra equations, TRS-ODEN and HODEN improve upon the baseline Neural ODE, but have slightly higher errors than our discovered equations with Lie symmetry regularization.\n\nWe agree that our work and TRS-ODEN cover different important aspects of ODE symmetries for different tasks. We will cite this work and explain the connections in the revised paper.\n\n> Effect of noise levels in Sec 5.1\n\n**Figure 4** in the supplementary pdf shows the equation discovery statistics under smaller noise levels. It is observed that the baseline without symmetry deteriorates quickly as the noise increases, while our symmetry-based methods constantly discover the correct equations.\n\n> If we don’t know the symmetry in the first two examples, what will be the result?\n\nIn this case, we need to discover the symmetry first. We use LieGAN [1] (without the autoencoder) to learn the linear symmetry in these examples. For the damped oscillator, we learn the symmetry $v=0.97x_2\\partial_1-x_1\\partial_2$. The resulting equivariant basis is close to the one corresponding to the true rotation generator. Solving this symmetry constraint gives similar equation discovery results (**Figure 4** left, red line). For the growth model, the learned symmetry is $v=1.95x_1\\partial_1+x_2\\partial_2$. Solving the constraint in eq (10) also results in 3 close-to-zero singular values and we obtain exactly the same equivariant basis as in Figure 2 (lower) in the paper. Therefore, the equation discovery results are also the same.\n\n[1] Generative Adversarial Symmetry Discovery. ICML 2023.\n\n> High-dimensional examples\n\nIn the supplementary PDF (**Table 4**), we provide two examples of learning higher-dimensional dynamics by solving the hard constraint (SEIR) and symmetry regularization (Lorenz).\n\nThe SEIR equations [1] model the evolution of pandemics by the number of susceptible, exposed, infected, and recovered individuals. It is a 4-dimensional ODE system with quadratic terms. We consider the following equations:\n\n$\n\\\\left\\\\{\n\\begin{aligned}\n\\dot S &= 0.15-0.6SI\\\\\\\\\n\\dot E &= 0.6SI-E\\\\\\\\\n\\dot I &= E-0.5I\\\\\\\\\n\\dot R &= -0.15+0.5I\n\\end{aligned}\n\\\\right.\n$\n\nWe assume an intuitive symmetry of scaling the number of recovered individuals proportional to the total population: $v=(S+E+I+R)\\partial_R$. We solve the constraint w.r.t this symmetry as in Sec 4.1. **Figure 2** (supplementary PDF) shows the reduced parameter space (from 60D to 34D).\n\nFor the 3D Lorenz system [2], we use parameters $\\sigma=0.5,\\beta=1.0,\\rho=0.0$. Similar to Section 5.3, we discover the symmetry first and use the learned symmetry to regularize equation discovery.\n\nPlease let us know if you have additional questions. If your concerns have been addressed, may we kindly request you increase the score for our paper?", "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": "19Zc3YFPkz", "reviewer": "Reviewer_zT4E", "summary": "The authors consider an estimation technique similar to the well known SINDy technique for recovering interpretable parameterizations of ordinary differential equations (ODEs). The authors primarily consider the context where there exists a Lie symmetry that constrains the solution space of ODEs to search over. The authors describe different ways in which to incorporate the symmetry information within the learning problem. The first is to consider linear constraints in which case the problem simplifies into another coefficient regression optimization. The second involves general constraints for which the authors propose a regularization scheme to constrain the function class. These methods are described in terms of both the case of observed data and a learned coordinate system. Finally, the authors consider a few different experiments where they try to estimate the coefficients for known data as well as predict the evolution of the data forward in time.", "questions": "Is there any existing theory on how well learned the symmetries can be from the data?\n\nHow does the computational cost of the method compare? I'm particularly interested in the increase in cost associated with the vector Jacobian product computation.\n\nIs there any intuition why the RMSE of W-SINDy is better than the proposed method in the case of successful identification?", "limitations": "The authors briefly discuss limitations regarding requiring known symmetries for the algorithm. It should probably be discussed more on how performance changes when using prior knowledge versus when learning these symmetries from the data.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Including the symmetry information is a natural approach to try to constrain the learning problem into something more manageable. The general SINDy learning problem can be difficult to interpret due to the requirement of the specification of the basis. In that sense, this method provides a good solution for imposing some constraints on the function space to make the problem more tractable.\n\nThe numerical results suggest the method performs very well compared to the baselines. For the forward prediction on the reaction diffusion system, the method also performs better than the related baselines. This shows the promise in some more general tasks such as time series forecasting  given sufficient regularity on the observed data.", "weaknesses": "The method largely focuses on cases where symmetries are linear. The regularization approach in equation 11 seems a bit arbitrary in its formulation (e.g. using the relative error) as well as it seems like the results are primarily applicable in the cases where the symmetry is linear.  \n\nSome of the empirical results are slightly confusing, in particular when the baselines fail catastrophically under some circumstances leading to high RMSE. However, in the cases where the the baselines successfully identify the equation, the baselines tend to have smaller RMSE. This may be a hyperparameter tuning issue, but it would seem more intuitive that the method with the symmetry should be easier to achieve a lower RMSE.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback. We address key comments below.\n\n> The method focuses on cases where symmetries are linear. The results are primarily applicable in the cases where the symmetry is linear.\n\n**This is a misunderstanding.** Our regularization method in Section 4.2 applies to general Lie symmetries, including linear and nonlinear ones. Specifically, the infinitesimal symmetry $v$ in eq. 11 can be either linear or nonlinear functions. Also, experimental results in Section 5.3 show that our method can work with nonlinear symmetry. Please let us know if you have additional questions regarding this point.\n\n> The relative error in eq. 11 seems arbitrary in its formulation\n\nUsing relative error in eq 11 is not arbitrary. We have explained the rationale for this formulation in the text following eq 11. The main idea is that an absolute loss would decrease with the magnitude of $h$. If we used the absolute error instead of the relative error, it would bias the discovered $h$ (towards a lower magnitude) and lead to a large coefficient estimation error.\n\nTo demonstrate this, we provide an ablation study on the Lotka-Volterra system, using both relative error and absolute error. From the table below, the absolute error causes a larger negative bias (computed as the mean of $(\\hat \\theta_i - \\theta_i) / \\theta_i$), meaning the discovered coefficients tend to have smaller magnitudes. It also increases the RMSE and reduces the success probability of finding the correct equation form. For reference, the results in the “None” and “Relative” rows correspond to Table 4, “L-V: SINDy” and “L-V: EquivSINDy-r” rows, in the paper.\n\n> Interpretation of empirical results, in particular when baselines achieve lower RMSE\n\nOur symmetry-based methods mainly increase the **success probability** of finding correct equation forms, which is the primary metric as we discussed in L275-277. This demonstrates the main benefits of using symmetries.\n\nSymmetry-based methods can have higher parameter estimation RMSE when the symmetry learned from data is not perfectly accurate. In comparison, when we know the exact symmetry in Section 5.1, the RMSEs of our method and SINDy w/o symmetry are similar.\n\nTo understand this, note that symmetry improves equation discovery by identifying a subspace of possible equations. When we explicitly solve the constraint (Sec 4.1), we enforce the model to search strictly within this space; when we use regularization (Sec 4.2), we encourage it to search around it. In either case, this smaller subspace makes it more likely to identify the correct solution. This explains the higher success probability of symmetry-informed methods. On the other hand, once the model (either with or without symmetry) manages to reach the proximity of the correct solution (i.e. it has found the correct equation form), symmetry no longer helps. At this stage, a slightly inaccurate symmetry learned from data may even cause a small bias, which explains why our symmetry-based method in Sec 5.2 sometimes has a higher RMSE. **Figure 1** in the supplementary PDF demonstrates the above argument with the equation space abstracted into a 2D plane.\n\nMoreover, to reduce the influence of inaccurate learned symmetry, we can refine the discovered equation w/o regularization. Specifically, we first perform the same experiment as in the paper with symmetry regularization. Then, we remove the regularization, fix the form of the equation, and fit the equation parameters under the equation fitting loss (same as baseline). The table shows that this refinement process can further reduce the parameter estimation error on the glycolytic oscillator. For reference, the “SINDy” and “EquivSINDy-r” rows here correspond to the “Sel’Kov” rows and RMSE (successful) columns in Table 4 of the paper.\n\nPlease let us know if the above discussions clear up the confusion. We will incorporate these comments in the revision.\n\n> Any existing theory on how well learned the symmetries can be?\n\nThat is an interesting question! To our knowledge, there isn’t such a theory that bounds the error of the learned symmetry in the latest works on symmetry discovery. Some relevant theoretical results include the necessary conditions for symmetry discovery with GAN [1], and quantitative metrics for learned symmetry [2,3]. We will leave this to our future work.\n\n[1] Generative Adversarial Symmetry Discovery. ICML 2023.\n\n[2] Latent Space Symmetry Discovery. ICML 2024.\n\n[3] LieGG: Studying Learned Lie Group Generators. NeurIPS 2022.\n\n> Computational cost of our method\n\nOur method moderately increases the computational cost. Excluding common procedures across methods, e.g. data loading, SINDy has an average run time of 11.62s and ours has an average of 20.44s over 50 runs on the L-V system, adding only <1x computational overhead. The computational cost also depends on how complicated the symmetry is. E.g. if we use an arbitrary simple symbolic function as symmetry (which certainly affects accuracy, but it’s fine since we are investigating efficiency), instead of the 5-layer MLPs in LaLiGAN, the average running time decreases to 13.07s, because the JVP computation is cheaper.\n\n> Any intuition why the RMSE of W-SINDy is better?\n\nWe have discussed this briefly in Appendix C.1 (L633-639). The intuition is that WSINDy computes loss based on the integration within a time interval, instead of time derivative errors at individual points. The white noise is averaged out over this time interval. As a discrete analogy, let $X_1,...,X_n$ be i.i.d standard Gaussian RVs, then their mean has reduced variance $1/n$. This explains why WSINDy is better at parameter estimation with our highly noisy data.", "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": "eFNrw5Mi6k", "reviewer": "Reviewer_HCHE", "summary": "This paper propose to leverage the symmetry to discover underlying dynamics (especially the ones described by autonomous ODEs) correctly from data. Specifically, the proposed method combines conventional symbolic regression approaches like SINDy with the symmetry-based constraints or regularizations. The used symmetry-based regularization is derived from the Lie group symmetry, namely the flow of the given dynamics should be $G$-equivariant with respect to the given Lie group $G$ of symmetries. For the case that the group of symmetry transformations is unknown, the authors use some recent techniques that can learn the unknown symmetry from data, and parameterize the Lie algebra of the learned symmetry with neural networks. In this case, the extracted symmetry regularization can be encoded by using the infinitesimal version of $G$-equivariance with the learned Lie algebra generators. The authors validate their proposed approach with some dynamical systems.", "questions": "The author mentions that though there are four forms of Lie point symmetry (i.e., equations (4 – 7) or equations (34 – 17)), they are empirically similar in terms of the performance. What would happen if the authors use multiple ones simultaneously instead of just one of these four? For example, how would the performance be if you use both (34) and (35) at the same time? I think this would work similarly to the Sobolev training [2], thus potentially improving performance, especially when the dataset is noisy.\n\n[2] Czarnecki et al. Sobolev Training for Neural Networks. NeurIPS 2017.", "limitations": "The authors mention that their method only consider the time-independent point symmeties for autonomous ODEs. They provides some potential future research directions for generalizing their approaches, e.g., for non-autonomous ODEs.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "**Significance.** Using symmetries to generalize the machine learning dynamics is very relevant topic. The combination of the symbolic regression with the symmetry regularization for recovering the underlying dynamics is thus very principle.\n\n***\n\n**Clarity.** The paper is well-written and easy to follow. The main paper is concise and highlights the advantages of the proposed technique without heavy math/details. Readers who want to know the detailed background and techniques on the proposed method can find the details in the appendix.", "weaknesses": "**Novelty.** The proposed technique is straightforward, if one already knows the group of symmetry transformations that should be applied. It is basically the regularization for the (infinitesimal) $G$-equivariance of the given Lie group $G$, i.e., $\\lVert f(g \\* x) - g \\* f(x) \\rVert$. Therefore, the interesting part of the proposed method is the case that the symmetry is unknown. In this case, the authors rely on [1], which can find the infinitesimal generators of Lie group of symmetries from data. Because they use [1] without clear modification/enhancements, there is some doubt as to whether the combination of two already known techniques (equivariance regularization and symmetry discovery) can be sufficiently novel.\n\n[1] Yang et al. Latent Space Symmetry Discovery. arXiv 2023.\n\n***\n\n**Experiments.** The proposed method is evaluated on relatively simple dynamics. It would be beneficial if the authors could provide experimental results on more complicated dynamics, such as chaotic systems.\n\n***\n\n\nThere are some typos, e.g., the denominator of (35), and the summation spaces ($v$ instead of $g$) of (35) and (36)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback and the recognition of our paper's significance and clarity. We address key comments below.\n\n> Novelty\n\nOur method is *not* a combination of existing techniques.\n\nFirst, we develop the pipeline for using different kinds of symmetries in equation discovery. For a known linear symmetry, we propose to solve the symmetry constraint explicitly with Lie’s infinitesimal criterion (Section 4.1), which is our novel contribution.\n\nThen, when the symmetry is unknown, we use symmetry discovery techniques to learn it from data. It does not have to be LaLiGAN, but we choose it in our experiments because of its ability to express nonlinear symmetries. We also make important adaptations to incorporate the discovered symmetry into our approach. E.g. We extract the infinitesimal action in the original state space as in eq (12), which is essential for symmetry regularization. This is not discussed in LaLiGAN (which only addressed linear infinitesimal actions in latent space); our approach opens up new possibilities for applying the learned symmetries. We also proposed the relative losses for regularization to prevent bias in the learned equations. An ablation study for relative loss is provided in the supplementary pdf (Table 1). We show that relative loss improves success probability and reduces parameter estimation errors on Lotka-Volterra equations.\n\n> Experiment with more complicated dynamics\n\nIn the supplementary PDF (**Table 4**), we provide two examples of learning higher-dimensional dynamics by solving the hard constraint (SEIR) and symmetry regularization (Lorenz).\n\nThe SEIR equations [1] model the evolution of pandemics by the number of susceptible, exposed, infected, and recovered individuals. It is a 4-dimensional ODE system with quadratic terms. We consider the following equations:\n\n$\n\\\\left\\\\{\n\\begin{aligned}\n\\dot S &= 0.15-0.6SI\\\\\\\\\n\\dot E &= 0.6SI-E\\\\\\\\\n\\dot I &= E-0.5I\\\\\\\\\n\\dot R &= -0.15+0.5I\n\\end{aligned}\n\\\\right.\n$\n\nWe assume an intuitive symmetry of scaling the number of recovered individuals proportional to the total population: $v=(S+E+I+R)\\partial_R$. We solve the constraint w.r.t this symmetry as in Sec 4.1. **Figure 2** (supplementary PDF) shows the reduced parameter space (from 60D to 34D).\n\nFor the 3D Lorenz system [2], we use parameters $\\sigma=0.5,\\beta=1.0,\\rho=0.0$, leading to non-chaotic dynamics. Similar to Section 5.3, we discover the symmetry first and use the learned symmetry to regularize equation discovery.\n\nThe presence of a nontrivial Lie symmetry often implies some form of simplification or integrability [3], which contradicts the nature of chaotic systems. We are not aware of any nontrivial Lie symmetry in chaotic dynamics.\n\n[1] Abou-Ismail A. Compartmental Models of the COVID-19 Pandemic for Physicians and Physician-Scientists. SN Compr Clin Med. 2020;2(7):852-858.\n\n[2] Brunton et al. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 2016\n\n[3] Olver. Applications of Lie groups to differential equations. Springer Science & Business Media, 1993.\n\n> Equations (34-37)\n\nThank you for noticing these! These are **not** typos but do require more clarification. The summation spaces should indeed be $v \\in B(\\mathfrak g)$, because we can only enumerate the finite list of Lie algebra basis components, not the infinite group elements. The group elements in (35)(36) are obtained from $v$ by $g=\\exp(\\epsilon v)$, where $\\epsilon$ is chosen manually. We can make the dependency of $g$ on $v$ explicit to prevent confusion.\n\nAlso, the denominator of (35) should indeed be the difference between $f(gx)$ and $f(x)$. If we use $f(gx)$ instead, the scale of the denominator would also depend on $x$. As a result, regions closer to the origin in the state space have larger weights (smaller denominators) during loss computation. It’s helpful to consider what would happen in the limit of $\\epsilon \\to 0, \\tau \\to 0$. The formulas inside the norms in the numerator and denominator in (35) both become $O(\\epsilon\\tau)$, while $f(gx) = f(\\exp(\\epsilon v)x) = f(x+O(\\epsilon)) = x+O(\\epsilon)+O(\\tau)+O(\\epsilon\\tau)$, which does not match the numerator.\n\n> Sobolev training\n\nThank you for the reference. The idea of Sobolev training is indeed relevant in this case. In the supplementary pdf (**Table 3**), we did additional experiments on both systems in Section 5.3, where we tried the combinations of (34)+(35) and (35)+(36). The latter involves $h$ and $f=\\mathrm{odeint}(h)$, which is more similar to the formulation of Sobolev training because we are learning the dynamics $h$ instead of the group action $v$. We compared them with the best results from using a single loss term. The combination (35)+(36) has a slightly higher success probability of finding the correct equation forms than the original results using single losses on the Lotka-Volterra system. For reference, the full original results are in Table 5 in the paper.", "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": "VB1sAVLt2X", "reviewer": "Reviewer_CnKF", "summary": "The authors developed an approach that allows exploitation of symmetry in common equation discovery algorithms and test their approach on dynamic systems with and without symmetry.\nFor this they make use of Lie groups.", "questions": "I want to highlight a major question before all others.\nWith respect to the results in Table 2:\n- If there is no symmetry in the system, shouldn't your approaches fallback to their respective base approaches, i.e. SINDY and GP? Why does your approach still outperform them so consistently, despite there being no exploitable symmetry in the ODEs?\n- Did you look into the results of the LaLiGAN? Did you recover the non-existence of symmetries and _still_ performed better? Or did it find some symmetry which your approach then used?\n\nOther questions:\n\n- 78/79: Thank you for highlighting works that already exploit symmetries. Can you provide a brief description of how you differ from Loiseau & Brunton 2018 and Guan et al. 2021?\n- I have a minor question for Prop 3.2. What does it mean for the flow map  $f_\\tau$ to be equivariant to the G-action on X? Does that mean that $f_\\tau$ needs to still be a valid flow map after an element of the symmetry group has been applied to $X$?\n- For Theorem 3.3 I have a few minor questions:\n\t- What is $\\epsilon$  in the $g = \\exp(\\epsilon v)$ term? \n\t- Why is $g$ defined this way? Would it not be possible to calculate the Jacobian otherwise?\n\t- Is it guaranteed that $g \\in G$? Does that only hold because $G$ is a symmetry group? Or generally?\n\t- I assume that $v$ is a matrix, please correct me if I am wrong, what does it mean to perform $\\exp(\\epsilon v)$ in that case? Is this the matrix exponential?\n- What is $L_v$ in Proposition 4.1? Is it a linear transformation? \n- 164/165: I don't understand what $c$ is and what its purpose is\n- 180: The set of polynomials doesn't include products $x_1 x_2$ like in the example for SINDY in line 163, correct?\n- 183: Is $W$ vectorized with respect to rows $\\{[w_{11} ... w_{1n}], ...\\}$  or columns $\\{[w_{11} ... w_{n1}], ...\\}$?\n- Where does the factorization in Equation (10) come from? Are there established methods that easily generate this? What do its components mean?", "limitations": "In my opinion the only limitation of this work lies in the experiments in 5.3, as I stated before.\nI hope the authors can either point out how their approach differs from vanilla SINDy/GP and explain why it is better or demonstrate, perhaps using a \"neutral\" operator instead of LaLiGAN, that it does fall back and that the difference is due to LaLiGAN.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- Despite Lie groups being a new topic for me, I was able to follow the well written high-level explanations in Sections 3.1 and 3.2\n- I like the approach that was used to exploit additional structure in the search\n- The results are very strong on the chosen systems, though I am still left sceptic of _how much_ the results are better in Table 1.\n- The experiments section is written concisely.", "weaknesses": "_Weaknesses_\n- Just looking at Figure 1, I find it hard to recognize at which points the dynamics are searched, i.e. where the equation discovery itself is happening. From the text prior to Figure 1, I expect the authors approach to find the symmetry first, constrain the search space and then perform equation discovery, but Figure 1 suggests otherwise.\n- Chapter 4 could benefit from examples (or a running example) how corresponding functions might look like, to facilitate understanding. In particular since $M_\\Theta$ can be computed in a symbolic fashion. Since the work is already tightly packed, such examples can also be mentioned in the appendix.\n\n\n_Minor complaints/typos_\n- A small arrow/text indicating the start \"Start\" in the flowchart on the left side of Figure 1 would make orientation slightly easier.\n- Section 4 was harder to read, compared to section 3\n- 261: With RK4 you mean Runge-Kutta 4 methods?\n- Figure 2 could be bigger to improve readability.\n\n\n\nTo summarize: I enjoyed this work and find it relevant to the community and a good contribution. However, the explanations in Chapter 4 were partly unclear to me and could benefit from an additional revision.\nThe most conflicting point for me are the results of Table 2, which I hope to have clarified during the discussion.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the valuable feedback. We address key comments below, starting with the major question about Tab. 2.\n> *If there is no symmetry*, why does our approach still outperform baselines?\n\nThere *are* symmetries in the tasks in Tab. 2. If there were no symmetry, our approach would not outperform baselines. The confusion might arise because we mentioned “these systems do not possess any *linear* symmetries” in L348-349. But they do have *non-linear* symmetries. For the Lotka-Volterra system, LaLiGAN learned a “distorted” rotation symmetry in their paper. We used LaLiGAN and obtained a similar symmetry, which we then used for regularization.\n\n**This should also address the concerns in the Limitations part**. Please let us know if you have other questions about this point.\n> Understanding Fig. 1\n\nYour understanding of the text is correct. In Fig. 1, equation discovery is performed in the last column (green background, titled “learning”). The confusion might have arisen because we titled the first column “dynamics”, where we actually meant the data collected from dynamical systems. We will change the block titles to avoid potential confusion. Thanks for pointing this out!\n> Example of computing $M_\\Theta$\n\nWe had an example in Appendix B.2 (L605-608). We further explain it here.\n\nDefine $\\Theta(x_1,x_2)=[1,x_1,x_2,x_1^2,x_1x_2,x_2^2]^T$, and $f(x_1,x_2) = [x_2^2,x_2]^T$. Each row of $M_\\Theta(f)$ corresponds to one component of $f$. For the first row, we look at $f_1(x_1,x_2)=x_2^2$. It is the linear combination of the elements in $\\Theta$ with coefficients $[0,0,0,0,0,1]$, which becomes the first row of $M_\\Theta(f)$. Similarly, the second row is $[0,0,1,0,0,0]$.\n> On our difference from related works (Loiseau & Brunton 2018 and Guan et al. 2021)\n\nOur work addresses a broad range of symmetries in dynamical systems; our pipeline is also applicable when we do not know the symmetry a priori. These related works, however, only considered a few instances of symmetries based on prior knowledge about specific systems. E.g. Loiseau (2018) only considered the time translation symmetry (from the energy-preserving principle) of a reduced-order modeling of Navier-Stokes equations. Guan (2021) considered a set of reflection and permutation symmetries of the proper orthogonal decomposition (POD) coefficients.\n> Other questions\n\nThe reviewer raised questions about specific details in the paper. We believe all related contents in our paper are mathematically correct, but agree some explanation would be helpful. We answer the questions below and will add necessary details to the updated version.\n* L261: yes, we mean the Runge Kutta 4 method by RK4.\n* Prop 3.2: $f_\\tau$ is equivariant to $G$-action on $X$, if $f_\\tau(g.x)=g.f_\\tau(x),\\ \\forall g,x$. This is the definition of equivariance. It is somewhat inaccurate to state that “$f_\\tau$ is still a valid flow map after a group element is applied”. The group element $g$ does not transform the function $f_\\tau$ itself, but its input and output $x, f_\\tau(x) \\in X$.\n* Thm 3.3: this involves some knowledge of Lie's theory. We refer to Appendix A.1 for background knowledge and answer the specific questions below.\n    * $\\epsilon\\in R$ is the scaling factor that determines the “size” of the transformation $g$.\n    * $g$ is defined as such because we consider the Lie group of continuous transformations. It is not related to the Jacobian computation. If we have a Lie algebra element $v$ (or $\\epsilon v$, which is still in the Lie algebra), we can always get a group element $g\\in G$ by applying the exponential map.\n    * (Is it guaranteed that $g\\in G$?) Yes, it follows from definitions of Lie group and Lie algebra. See above.\n    * (Does $g\\in G$ hold because $G$ is a symmetry group?) No, it holds for a general Lie group. This is how we describe a Lie group element: by its counterpart in the Lie algebra, a flat vector space that is easier to deal with. At this point, we don’t need to consider whether $G$ is the symmetry of ODE. We just consider $G$ as a Lie group itself.\n    * (Is $v$ a matrix? What is $\\exp$?) $v$ could be a matrix, but not necessarily. More generally, $v$ could be a nonlinear vector field, and $\\exp$ can be understood as following the flow generated by that vector field. Intuitively, think of a 2D plane. The vector field $v: R^2 \\to R^2$ consists of arrows pointing in various directions at every point. Exponentiating the vector field is like starting at one point and walking in the arrow directions for a certain time $\\epsilon$. The path you take following these arrows traces your journey, and the ending point represents the effect of the vector field exponentiation. In terms of computation, if $v$ is a general vector field, $\\exp$ is computed by integrating the vector field for a certain time interval; if $v$ is a matrix, it is equivalent to matrix exponential.\n* Prop 4.1: yes, $L_v$ is a matrix, and $L_v x$ is mat-vec multiplication.\n* L164-165 (what is $c$?): it reads as: if a function $f$ is a linear combination of the $p$ functions in $\\Theta$, with the coefficients $c\\in R^p$, i.e. $f(x)=c^T\\Theta(x)$, then $M_\\Theta(f)$ is defined as $c$.\n* L180 (does the set of polynomials include $x_1x_2$?): Yes, it includes products of different $x_i$’s. Polynomials of multiple variables include terms like $x_1^{q_1}...x_d^{q_d}$, where the polynomial degree is defined as $q = \\sum_i q_i$.\n* L183: As in common practice, matrix vectorization stacks the columns: $\\mathrm{vec}(W)=[w_{11},...,w_{d1},w_{12},...]$.\n* Eq (10) refers to singular value decomposition. There are built-in SVD functions in any modern numerical computation package (e.g. torch.svd). Of particular interest is $Q^T$, the right singular vectors corresponding to zero singular values. They span the null space of $C$, i.e. the solution space of $W$.\n\nPlease let us know if you have other questions. If your concerns have been addressed, may we kindly request you increase the score?", "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:
We thank the reviewers for their detailed and valuable feedback. We are encouraged that they find our work to be a clearly motivated idea (R2,3) towards the interesting problem of equation discovery (R1). We are also glad that they find our paper well-written and easy to follow (R1,4), our method principled and theoretically grounded (R2,4), and demonstrated to be effective in a wide range of tasks (R2,3).

## Common questions

Here we answer some common questions from the reviewers. This is a high-level summary; more detail can be found in individual responses.

> On the connection and difference from LaLiGAN

Our work has a completely **different goal** from LaLiGAN. Our method aims to discover equations using symmetry as an inductive bias. LaLiGAN aims to discover unknown symmetry. Only when we do not know the symmetry a priori, we use LaLiGAN as a tool to discover the symmetry first and use the discovered symmetry to regularize equation learning.

Also, our work has made important adaptations to incorporate learned symmetry into equation discovery, e.g. extracting the infinitesimal action in eq (12) and using relative regularization loss in eq (11).

> The types of symmetries in different experiments

There is some confusion about what kind of symmetry we are using for each experiment. In particular, in Sec 5.3, we use *nonlinear* symmetries *discovered* by LaLiGAN for symmetry regularization. It is *not* true that the symmetry is linear (R2), or there is no symmetry in this case (R5).

We believe the above information is available in the paper. However, we will provide more clarifications in the revision to further avoid confusion.

## The supplementary PDF

Here we also provide an index of the supplementary PDF. It includes additional illustrations and tables of experiment results requested by the reviewers.
* Table 1: an ablation study for the relative regularization on the Lotka-Volterra equations. It is shown that the relative loss leads to a smaller parameter estimation error than the absolute loss.
* Table 2: refining symmetry-regularized solutions. The *learned* symmetries may introduce slight biases. This problem can be addressed by removing the symmetry regularization, fixing the equation form, and optimizing the equation parameters under the equation $L_2$ loss only. This process reduces the parameter estimation error.
* Table 3: using combined regularization terms can sometimes lead to better discovery. This is motivated by the idea of Sobolev training, thanks to the reference from R4.
* Table 4: equation discovery on high-dimensional systems, as requested by R2 and R4.
* Figure 1: an intuitive demonstration of how symmetry helps the search in the equation parameter space. The equation parameter space is abstracted into a 2D plane, where the symmetry constraint leads to a smaller 1D subspace and makes optimization easier.
* Figure 2: the equivariant basis for the SEIR model in Table 4. In this case, symmetry reduces the 60D parameter space to 34D, significantly reducing the complexity of the optimization problem.
* Figure 3: experiments in Sec 5.1 with different noise levels.
* Figure 4: Comparison with Neural-ODE-based methods, including TRS-ODEN, as requested by R2.

Meta-review: This paper enhances equation discovery by incorporating symmetry constraints from Lie groups, improving accuracy and simplicity. The method, tested on various dynamical systems, shows superior performance over baseline approaches by leveraging symmetry to refine and constrain the search space. Some confusion about the difference to LaLiGAN (which also proposes a simpler equation discovery method) can and should be clarified in the rebuttal and should be made very clear in the camera ready version. The rebuttals have resolved other reviewer concerns, and the proposed changes are straightforward for the final version. I concur with the reviewers' consensus on acceptance. Ensure the camera-ready version also clearly and early explains that the linear case restricts the solution space while the non-linear case introduces an inductive bias via symmetry regularization, as discussed with Reviewer zT4E.
Decision: Accept (poster)