Symmetry-Informed Governing Equation Discovery

Authors: jianke  yang
wang  rao
nima  dehmamy
robin  walters
rose yu ucsd
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 21:58:46

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.

Main paper: 
 introduction: 
 Discovering governing equations that can describe and predict observational data from nature is a primary goal of science. Compared to the black-box function approximators common in machine learning, symbolic equations can reveal a deeper understanding of underlying physical processes. In this paper, we consider a first-order dynamical system governed by an autonomous ordinary differential equation (ODE):EQUATIONwhere t ∈ T denotes time, x(t) ∈ X is the system state at time t, X ⊆ R d is the phase space of the ODE, ẋ(t) the time derivative, and h : X → R d 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 a 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 to 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 to 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 a 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 unnecessarily large space of equations and discover equations that are overly complex or do not conform to 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 to design equivariant networks (Bronstein et al., 2021; Finzi et al., 2021; Wang et al., 2021) , augment training data (Benton et al., 2020) , or regularize the model (Akhound-Sadegh et al., 2023; Otto et al., 2023) , leading to improved prediction accuracy, sample efficiency, and generalization.Specifically, we derive the equivariance constraints from the symmetries of ODEs and solve the constraints explicitly to compress the search space of equations. Alternatively, we can also use symmetry as a 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, we also incorporate a method for learning symmetries with nonlinear group actions (Yang et al., 2023a) . 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 time-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 a symmetry regularization term when the symmetry constraint cannot be explicitly solved. • In experiments across many dynamical systems with substantial noise, our symmetryinformed approach achieves higher success rates in recovering the governing equations and better robustness against noise.2 Related Works 
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 (Cranmer et al., , 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 et al., 2016) , these methods assume that the equation can be written as a linear combination of some predefined functions and optimize the coefficients with sparsity regularization. Subsequently, Champion et al. (2019) combines the sparse regression technique with an autoencoder network to simultaneously discover coordinate systems and governing equations. Weak SINDy (Messenger & Bortz, 2021) 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 as 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) .Symmetry plays a crucial 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 Finzi et al., , 2021;; Bronstein et al., 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 learning in dynamical systems governed by differential equations through data augmentation (Brandstetter et al., 2022) , symmetry loss (Huh et al., 2020; Akhound-Sadegh 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., 2021) , rotations and translations in space (Ridderbusch et 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 a priori. Some recent works have aimed to automatically discover symmetries from data (Liu & Tegmark, 2022; Yang et al., 2023b,a; Otto et al., 2023; van der Ouderaa et al., 2024) . When there is no available prior knowledge about symmetry, We use the adversarial framework in Yang et al. (2023a) to learn unknown symmetries from data and discover the equations subsequently using the learned symmetry. 
symmetry of odes: 
 Our goal is to discover the governing equations of a dynamical system (1) in symbolic form. Our dataset consists of observed trajectories of the system, {x(i) 0:T } N i=1 , where x (i)t ∈ X ⊆ R d denotes the observation at timestep t in the ith trajectory. We assume the observations are collected at a regular time interval ∆t. We omit the trajectory index i and consider one trajectory, x 0: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 a time-independent 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.The Lie point symmetry of an ODE is a transformation that maps one of its solutions to another. The symmetry transformations form the symmetry group of the ODE. In this paper, we focus on symmetry transformations that act solely on the phase space X ∋ x, without changing the independent variable t. We refer to these as time-independent symmetries, formally defined as follows.Definition 3.1. A time-independent symmetry group of the ODE (1) is a group G acting on X such that whenever x = x(t) is a solution of (1), x = g • x(t) is also a 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 to derive equivariance constraints on equations.Given a fixed time τ , the flow map associated with the ODE (1) is denoted f τ : X → X, which maps a starting point in a trajectory x(t 0 ) to an endpoint x(t 0 + τ ) after moving along the vector field h for time τ . Formally, f τ (x 0 ) = x(τ ), where x(t) is the solution of the initial value problemẋ(t) = h(x), x(0) = x 0 , t ∈ [0, τ ].(2) Proposition 3.2. Let G be a group that acts on the phase space X of the ODE (1). G is a symmetry group of the ODE (1) in terms of Definition 3.1 if and only if for any τ ∈ T , the flow map f τ is equivariant to the G-action on X.As our dataset consists of trajectories measured at discrete timesteps, consider the discretized flow map x t+1 = f (x t ), where f := f ∆t depends on the sampling step size ∆t of the data. From Proposition 3.2, this function of moving one step forward in the trajectory is G-equivariant, i.e.x t+1 = f (x t ) ⇒ g • x t+1 = f (g • x t ).(3)We can use the equivariance condition (3) to 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-independent symmetry group of the ODE (1). Let v ∈ g be an element in the Lie algebra g of the group G. Consider the flow map f defined in (2) for a fixed time interval. Denote J f and J g as the Jacobian of f (x) and g • x w.r.t x. For all x ∈ X and g = exp(ϵv) ∈ G, the following equalities hold:EQUATION)To understand Theorem 3.3, (4) 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 is infinitesimally transformed by v, and v(f (x))indicates the amount of change caused by the infinitesimal action on the output of f .2 The difference between these two terms vanishes if the function is equivariant. On the other hand, as both h and v 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(ϵv). Finally, we can 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 B.1.Theorem 3.3 provides the theoretic basis for our symmetry-based equation discovery methods in Section 4. We can enforce the equality constraints in Theorem 3.3 by solving the constraints explicitly (Section 4.1) or adding them as regularization terms (Section 4.2). 
symmetry-informed governing equation discovery: 
 Next, we discuss how to 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 1 . We introduce the components in the following subsections.We start from the case of linear symmetries. Many ODEs exhibit symmetries with linear actions on the state variable x. 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 of basis functions as in SINDy sparse regression (Brunton et al., 2016) , the infinitesimal formula (7) becomes a 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 Θ(x) ∈ R p as:EQUATIONwhere W ∈ R d×p is the learnable coefficient matrix. In SINDy, the function library Θ(x) is predefined. For example, Θ can be the set of all polynomials up to second order for a 2-dimensional system, i.e. Θ(x 1 , x 2 ) = [1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 ]T . Our goal is to solve for the parameter W that makes (7) hold for any x. To this end, we define the symbolic map M Θ : (R d → R) → R p , which maps a function to its coordinate in the function space spanned by Θ. Specifically, M Θ (f j ) = c for f j (x) = c T Θ(x), c ∈ R p . For a multivariate function f : R d → R n , we compute M Θ for each of its components and stack them into M Θ (f ) ∈ R n×p . A concrete example is provided in Appendix B.2. We have the following proposition. Proposition 4.1. Let G be a time-independent symmetry group of the ODE ẋ = W Θ(x) with linear actions. Let v ∈ g be an Lie algebra element with action X by v(x) = L v x, L v ∈ R d×d . Then, L v W = W M Θ (J Θ (•)L v (•)), where J Θ denotes the Jacobian of Θ. ℒ symm = ∥[v, ĥ ]∥ 2 Equation learning ⊙ [ ] Q 1 Q 2 [ ] β [ ] ĥ (x) = WΘ(x) = Θ(x) Symmetry v ∈ 𝔤 a b c a b c Data x (i) 0:T Enc. EquivSINDy L v = [ 0 1 -1 0] v = x 2 ∂ 1 -x 1 ∂ 2 arg min β ∥ • x -ĥ (x)∥ 2 Discover symmetry W Q Q : Equivariant basis β : Learnable parameters v vĥ ĥ• x = ĥ (x)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 or symmetry discovery techniques. We then enforce the symmetry by solving a set of constraints when possible and otherwise promote the symmetry through regularization.Proof. Substituting (8) into (7), we obtainJ v (x)W Θ(x) = W J Θ (x)v(x), ∀x ∈ X. As v acts on X linearly, i.e. v(x) = L v x, the constraint is equivalent to L v W Θ(x) = W J Θ (x)L v x.We then apply M Θ to both sides of the equation. Obviously, M Θ (Θ) = I p . Since (7) is true for all x, we haveL v W = W M Θ (J Θ (•)L v (•)).To calculate M Θ , we need to ensure that J Θ (x)L v x is still within the span of the functions in the given Θ(x). This is true when Θ(x) is the set of all polynomials up to a certain degree, which we prove in Appendix B.2. Proposition 4.2. The components of J Θ (x)L v x ∈ R p can be written as linear combinations of the terms in Θ(x) if Θ(x) is the set of all polynomials up to degree q ∈ Z + .Once we have the infinitesimal action L v and an appropriate function library Θ, M Θ can be computed in a purely symbolic way using sympy (Meurer et al., 2017) . Denoting M Θ,Lv := M Θ (J Θ (•)L v (•)) and vectorizing the matrix W , we can rewrite the equation in Proposition 4.1 as the following linear constraint on W :(-M T Θ,Lv ⊕L v )vec(W ) = 0 (9) where ⊕ is the Kronecker sum: A ⊕B = A ⊗ I + I ⊗ B. We concatenate the constraints for all representations L vi of the Lie algebra basis {v i } c i=1 into a single matrix C and solve the constraint using the singular value decomposition:EQUATION)The coefficient matrix resides in the r-dimensional nullspace of C and can be parametrized as vec(W ) = Qβ, β ∈ R r . This equivariant parametrization significantly reduces the number of free parameters from d × p to r, leading to parsimonious equations and easier training in a compressed parameter space. We name this approach as EquivSINDy-c, where c stands for constraint.The above procedure uses linear symmetries to inform sparse regression. However, when the data is high-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 SymmetryThe above procedure, EquivSINDy-c, applies to linear symmetry with a proper choice of function library Θ. Solving the equality constraint in Proposition 4.1 for arbitrary group action is generally more challenging. Technically, when the infinitesimal action v is expressed in a closed form, we can still apply M Θ to both sides of the equation and obtainM Θ (v(W Θ(•))) = W M Θ (J Θ (•)v(•)).However, the function library Θ needs to be carefully chosen based on the specific action v to ensure that M Θ is evaluated on functions within the span of Θ(x). 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 M Θ is intractable.We propose an alternative approach that is universally applicable to 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 to the equation loss from the base equation discovery algorithm, e.g. L 2 error between the estimated time derivative and the prediction from the learned equation.We consider the following relative loss based on infinitesimal group action v and a finite-time flow map f τ of the equation as in ( 5):EQUATIONwhere B(g) denotes the basis of the Lie algebra, f τ is obtained by solving the initial value problem (2), and J fτ denotes the Jacobian of f τ .We use a relative loss because the scale of both terms in the numerator and their difference is subject to the specific f τ . In the extreme case when h = 0, both terms are zero because f τ (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, it is a 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), r standing for regularization.The knowledge of ODE symmetries is 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 et al., 2023b,a; Otto et al., 2023) 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 ≤ GL(k; R) as exp(ϵv)x = (ψ • exp(ϵL v ) • ϕ)(x), where the networks ψ and ϕ define an autoencoder and L v ∈ R k×k is the matrix representation of v ∈ g. Our dataset consists of ODE trajectories sampled at a fixed rate ∆t, each given by x 0:T . We can learn the equivariance of f : x t → x t+1 by feeding the input-output pairs (x i , x i+1 ) into LaLiGAN. We extract the infinitesimal action asEQUATIONwhere z = ϕ(x) is the latent mapping via encoder network, and the Jacobian-vector product can be obtained by automatic differentiation. This infinitesimal action is then used to compute the regularization ( 11).Alternatively, we can also discover equations in the latent space of LaLiGAN by solving the equivariance constraints in Section 4.1, as the latent dynamics z t → z t+1 has linear symmetry L v . 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 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. 
experiments: 
 The experiment section is organized as follows. In Section 5.1, we will demonstrate how to solve the linear symmetry constraints (Section 4.1) on some synthetic equations. In Section 5.2, we will apply the same equivariant model to discover equations in a symmetric latent space for some highdimensional observations. In Section 5.3, we consider two well-studied benchmark problems for dynamical system inference under significant noise levels, where we use the regularization technique (Section 4.2) based on learned symmetries.Data generation. For each ODE system, we sample N random initial conditions x 0 from a uniform distribution on a specified subset of X ⊆ R d . Starting from each initial condition, we integrate the ODE using the 4th-order Runge-Kutta (RK4) method and sample with a regular step size ∆t to get the discrete trajectory x 0:T . We use different N , T and ∆t for each system, as specified in Appendix D.Then, unless otherwise specified, we add white noise to each dimension of the state x at a noise level σ R . The scale of the noise depends on the variance in the data within each state dimension x i :σ i = σ R • std(x i ).We will report the noise level σ R for each system. After adding the noise, we apply a Gaussian process-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θ fi f i (x), where θ fi is a nonzero constant parameter and f i is a function of x. Also, we expand the discovered equation as θfi fi (x), where θfi ̸ = 0. The discovery is considered successful if all the terms match, i.e. {f i } = { fi }. The probability is computed as the proportion of successful runs. As each dynamical system consists of multiple equations (one for each dimension), we evaluate the success probability of discovering each individual equation (e.g. "Eq. 1", "Eq. 2" columns in Table 1 ) and the joint probability of successfully discovering all the equations (e.g. "All" column in Table 1 ). 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.The following metrics are also considered:• RMSE for parameter estimation ∥θ -θ(k) ∥ 2 /K, where K is the number of runs, θ = (θ f1 , θ f2 , ...) is the constant parameters in true equation terms as defined above, and θ(k) the estimated parameters of the corresponding terms from the kth run. If f i is not in the learned equation, θfi = 0. 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 set to 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 (Brunton et al., 2016) and Genetic Programming (GP) with PySR package (Cranmer, 2023) , and their respective variants using the variational formulation, i.e. Weak SINDy (WSINDy) (Messenger & Bortz, 2021) and D-CODE (Qian et al., 2022) . Based on different experiment setups and the guidelines in Figure 1 , we apply our methods referred to as EquivSINDy / EquivGP, where the suffix c stands for solving the linear symmetry constraint, and r stands for using the symmetry regularization.1 x1 x2 x 2 1 x1x2 x 2 2 h1(x) h2(x) 1 x1 x2 x 2 1 x1x2 x 2 2 h1(x) h2(x) 0.4 0.2 0.0 0.2 0.4 Damped Oscillator: Q 1 , Q 2 1 x1 x2 x 2 1 x1x2 x 2 2 h1(x) h2(x) Growth: Q 1, 2, 3Figure 2 : Solutions for equivariant constraints of (13, 14). In ( 13), the 2D parameter space is spanned by Q 1 and Q 2 . In ( 14), the 3D parameter space is spanned by Q 1,2,3 , each marked with a different color.In this section, we consider two dynamical systems which possess some linear symmetries. These symmetries are commonly observed across many systems and their presence can be easily detected. Therefore, we assume we know these symmetries and focus on symmetry-informed equation discovery by solving the symmetry constraints.Damped Oscillator. This is a two-dimensional system with rotation symmetry, characterized by the in-EQUATIONGrowth. This system exhibits a scaling symmetry under(x 1 , x 2 ) → (a 2 x 1 , ax 2 ). The correspond- ing infinitesimal generator is v = 2x 1 ∂ 1 + x 2 ∂ 2 .It is governed by ( 14). EquivSINDy-c 1.00 1.00 1.00Figure 2 visualizes the parameter spaces for these two systems under the equivariance constraint. The original parameter space is W ∈ R 2×6 as we build the function library Θ with up to second-order polynomials. For (13), the equivariance constraint reduces the parameter space to two dimensions. The equivariant basisQ 1 = Qe 1 and Q 2 = Qe 2 are visualized.For ( 14), the parameter space is reduced to three dimensions. Each equivariant basis component only affects one equation term, visualized in three distinct colors. We can observe that the equivariant models have much fewer free parameters than non-equivariant ones.Table 1 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.1. Overall, these results fully demonstrate the advantage of applying symmetry to compress the parameter space.Figure 3 : Reaction-diffusion system and the prediction error in the high-dimensional input space using equations from different methods. The means and standard deviations (shaded area) of errors over 3 random runs are plotted. We also show an example of discovering equations in the latent space with an equivariant model. We consider a λ -ω reaction-diffusion system visualized in Figure 3 . The system is sampled at a 100 × 100 spatial grid, yielding observations in R 10000 . See Appendix D.1 for more details. We apply LaLiGAN (Yang et al., 2023a) to discover a latent space R 2 where the dynamics z t → z t+1 is equivariant to a linear group action. Meanwhile, we use equivariant SINDy with the learned symmetry to discover the equation for the latent dynamics. We compare our approach with SINDy Autoencoder (Champion et al., 2019) without any symmetry and LaLiGAN + SINDy, i.e. learning the equations in the LaLiGAN latent space but without equivariance constraint.Figure 3 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. 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-Volterra System describes the interaction between a predator and a 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.EQUATIONGlycolytic Oscillator (Sel'Kov, 1968) is a biochemical system of two ODEs with complex interactions in cubic terms. We use the same constants in the equations as in Qian et al. (2022) .ẋ1 = 0.75 -0.1x 1 -x 1 x 2 2 ẋ2 = 0.1x 1 -x 2 + x 1 x 2 2 (16)Since we assume no knowledge of the symmetries here, we first run LaLiGAN (Yang et al., 2023a) to 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 errors averaging over all random experiments.In Appendix C.2, we provide additional results to show that the equations discovered by our method also achieve lower prediction error. We also compare the different regularization options based on Theorem 3.3. 
conclusion: 
 We propose to incorporate symmetries into equation discovery algorithms. Depending on whether the action of symmetry is linear, we develop a pipeline (Figure 1 ) with various techniques, e.g. solving the equivariance constraint and using different forms of symmetry regularization, to perform symmetry-informed 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 to 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, our method can be generalized to other settings, such as time-dependent symmetries and non-autonomous ODEs. We provide a more general problem definition and discuss these potential generalizations in Appendix A.2. However, in these problems, it is less likely that we know the symmetries before discovering the equations. Extending our symmetry-informed methodology to these scenarios would require further breakthroughs in learning symmetries for general differential equation systems, which could be an interesting future research direction. 
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] Justification: 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, we report the mean and standard deviation over all the replicated runs.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, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are 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).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [No] Justification: We did not keep track of the execution time as most experiments are solving low-dimensional problems and can be done within minutes on a single CPU or GPU.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 9. 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 no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The data used in the paper are generated by numerical methods and have no safety risks.Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS 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. 
 
 Appendix: We use the theory of Lie groups to describe the symmetry transformations of ODE solutions. This section provides a brief overview of some fundamental concepts of Lie group and its application to differential equations.A.1 Lie Groups Lie groups. A Lie group is a smooth manifold with a compatible group structure. We use Lie groups to describe continuous transformations. For example, all n × n invertible matrices with real entries form the general linear group GL(n; R) with matrix multiplication as the group operation.The Lie algebra g := T e G is defined as the tangent space at the identity element e of the group G.For example, The Lie algebra of general linear group GL(n, R) consists of all real-valued matrices of size n × n. We can use a vector space basis {v i ∈ g} to describe any Lie algebra element as v = ϵ i v i , where ϵ i ∈ R. Lie algebra can be thought of as the space of infinitesimal transformations 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 a map.A group action defines how abstract group elements transform other objects in a space V . The group action is a function Γ : G × V → V which maps the identity e to identity transformation, i.e. Γ(e, v) = v, and is compatible with group element composition, i.e. Γ(g 1 , Γ(g 2 , v)) = Γ(g 1 g 2 , v). We can also think of the group action as mapping a group element to a function V → V and write Γ g := Γ(g, •). The group action Γ induces an infinitesimal action γ of its Lie algebra given by γv (x) = d dϵ (Γ exp(ϵv) (x))| ϵ=0 .For a smooth manifold V , the infinitesimal action can be interpreted as a vector field on V that generates the flow coinciding with the group action of exp(ϵv). Given the coordinate (x 1 , ..., x d ), the vector field can be written asγ v = i v i ∂ i , where ∂ i ≡ ∂∂xi and v i : V → R are the components of γ v . A group action is called a representation if V is a vector space and Γ g is a linear function on V . When V = R n , the representation Γ g ∈ GL(n; R) is an invertible n × n matrix that transforms v ∈ V by matrix multiplication. Correspondingly, the induced Lie algebra representation γ v ∈ gl(n, R) is also an n × n matrix. In the following discussion, we will refer to 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, we use g ∈ G to refer to both the group element and the group action Γ g and write g • x := Γ g (x). Similarly, we use v ∈ g to refer to both the Lie algebra element and its infinitesimal action and write v(x) := γ v (x).In this work, we are mainly concerned with the autonomous ODEs in the form (1). However, to introduce the general notion of Lie point symmetry of ODEs, it is helpful to consider a more generalized setup as follows:EQUATIONwhere t ∈ T ⊆ R denotes time, i.e. the only independent variable in ODEs, and x ∈ X ⊆ R d denotes the dependent variables.To formalize the concept of symmetries, we take a geometric perspective on the differential equations (17). The independent and the dependent variables form a total space E = T × X. We define the firstorder jet space of the total space as J = T × X × X 1 . Its coordinates, (t, x 1 , ..., x n , x ′ 1 , ..., x ′ n ) ∈ 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.Then, we can represent the ODE (17) through the map ∆ : J → R n with componentsEQUATIONA symmetry of the system ∆ transforms one of its solutions to another. More specifically, we need to define the group action on a point in the total space E = T × X and its induced action on a function T → X.We express the action of group element g in the total space asEQUATIONwhere t and x are functions over T × X.To define the induced group action on the function f : T → X, we first identify f with its graphγ f = {(t, f (t)) : t ∈ Ω T } ⊂ T × X (20)where Ω T ⊂ T is the time domain. The graph of the function can be transformed by g asEQUATIONThe set gγ f may not be the graph of another single-valued function. However, by choosing a suitable domain Ω T , we ensure that the group elements near the identity transform the original graph into a graph of some other single-valued function x = f ( t) with the transformed graph γ f = g • γ f . In other words, the transformed function f = g • f is defined by the transformation of the function graph.Definition A.1 (Def 2.23, Olver (1993) ). A symmetry group of ∆ (18) is a local group of transformations G acting on an open subset of the total space T × X that whenever x = f (t) is a solution of ∆, and wheneverg • f is defined for g ∈ G, then x = (g • f )(t)is also a solution of the system.Let v be a vector field on the total space T × X with corresponding one-parameter subgroup exp(ϵv).The vector field describes the infinitesimal transformations to the independent and the dependent variables. It can be written asEQUATION)The vector field also has its induced infinitesimal action on the jet space T × X × X 1 described by the prolonged vector field:EQUATION)where ϕ(1)i can be computed by the prolongation formula (Theorem 2.36, Olver (1993)) as ϕEQUATIONWe denote ϕ(t, x) := [ϕ 1 (t, x), ..., ϕ d (t, x)] T as the vector field components in X.As a special case, for the time-independent symmetries we considered previously, ξ = 0 and ϕ(t, x) = ϕ(x), so we haveϕ (1) i = j ∂ϕ i ∂x j x ′ j (26)Another special case is time translation. When ξ = 1 and ϕ = 0, we have ϕ(1) i = 0.The following infinitesimal criterion gives the symmetry condition in terms of the vector field (and its prolongation).Theorem A.2 (The infinitesimal criterion. Thm 2.31, Olver (1993) ). Let ∆ be a system of ODEs of maximal rank defined over M ⊂ T × X, X ⊆ R d . A local group of transformations G acting on M is a symmetry group of the system ∆ iff3 for every infinitesimal genrator v of G, v (1) [∆ i ] = 0, i = 1, ..., d, whenever ∆ = 0.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 toEQUATIONFor a more comprehensive understanding of the application of Lie groups in differential equations, we refer the reader to Olver (1993) .Proposition B.1. Let G be a group that acts on X, the phase space of the ODE (1). G is a symmetry group of the ODE (1) in terms of Definition 3.1 if and only if for any τ ∈ T , the flow map f τ is equivariant to the G-action on X.Proof. The proof follows from applying various definitions. If G is a symmetry group of (1), then for any solution x = x(t), x = g • x is still a solution to (1). By definition of the flow map f τ , we have f τ ( x0 ) = xτ , where x0 := g • x(t; x 0 ) and xτ := g • x(t + τ ), for any x 0 ∈ X. That is, f τ (g • x 0 ) = g • x(t + τ ) = g • f τ (x 0 ), i.e. f τ is G-equivariant.On the other hand, for any solution to the ODE (1) x = x(t) where x(0) = x 0 , we consider the transformed function x = g • x. Because f τ is G-equivariant, we haveEQUATIONTaking the derivative w.r.t τ at τ = 0, we have d dt x(t) = h( x(t)), which means x is also a solution to the ODE (1).Theorem B.2. Let G be a time-independent symmetry group of the ODE (1). Let v ∈ g be an element in the Lie algebra g of the group G. Consider the flow map f defined in (2) for a fixed time interval. For all x ∈ X and g = exp(ϵv) ∈ G, the following equations hold:EQUATION)Proof. From Proposition B.1, the flow map f is equivariant to G, so (29) follows directly.Consider a group element in the identity component that can be written as g = exp(ϵv), v ∈ g.Taking the derivative with respect to ϵ on both sides of (29) at ϵ = 0, we haved dϵ f (g • x) ϵ=0 = d dϵ g • f (x) ϵ=0 ⇒ d dϵ f (g • x) ϵ=0 = v(f (x)) ⇒ J f (exp(0v)(x))v(x) = v(f (x)) ⇒ J f (x)v(x) = v(f (x))(31) and ( 32) can be proved similarly. Let f τ (x 0 ) = t0+τ t0 h(x)dt (where x(t 0 ) = x 0 ) be the flow map given the time interval τ . Note that d dτ | τ =0 f τ (x) = h(x). Since ( 29) is true for any τ , we can take the derivative with respect to τ on (29) at τ = 0 and obtain (31).Finally, we take the derivative with respect to τ on (30) at τ = 0. We haved dτ J fτ (x) τ =0 ij = ∂ ∂x j d dτ τ =0 f i τ (x) = [J h (x)] ij and d dτ v(f τ (x)) τ =0 = J v (f τ =0 (x))h(x) = J v (x)h(x)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 toEQUATIONor, in its matrix form, 30) is similar to some results in the existing literature, e.g. Section 4.2 in Finzi et al.J v (x)h(x) = J h (x)v(x), where v(x) = [ϕ 1 (x), ..., ϕ d (x)] T . Remark B.4. ((2021) and Theorem 3 in Otto et al. ( 2023). ( 30) is a generalization of these results which does not require the linearity of the function f or the infinitesimal action v.In Section 4.1, we have defined the symbolic map M Θ : (R d → R n ) → R n×p , which maps each component of a function with multivariate output to its coordinate in the function space spanned by Θ. This map is introduced to account for the fact that the constraint Proposition 4.1 holds for any x.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 = 2 and Θ the set of all polynomials up to second order, i.e. Θ(x 1 , x 2 ) = [1, x 1 , x 2 , x 2 1 , x 1 x 2 , x 2 2 ] T . Define f 1 (x 1 , x 2 ) = [x 2 2 , x 2 ] T and f 2 (x 1 , x 2 ) = x 3 1 .We have M Θ (f 1 ) = 0 0 0 0 0 1 0 0 1 0 0 0 , and M Θ (f 2 ) undefined because the cubic term cannot be expressed in terms of a linear combination of functions in Θ.Thus, for a given infinitesimal action v, we need to ensure that an appropriate function library Θ is chosen for SINDy so that it can express the equivariant functions to the given action. For a linear action v(x) = L v x, we prove that the polynomial functions satisfy the requirement.Proposition B.5. Let Θ(x) be the set of all polynomial functions up to degree q (q ∈ Z + ). Then, the components of J Θ (x)L v x ∈ R p can be written as linear combinations of the terms in Θ(x), so thatM Θ (J Θ (•)L v (•)) is well-defined.Proof. Consider the function library Θ(x) that includes all polynomial terms up to degree p ≥ 1. The Jacobian J Θ (x) captures the partial derivatives of each component of Θ(x), which are polynomials up to degree q -1. The linear transformation L v maps x to L v x, and therefore, is a first-degree polynomial in terms of x. Hence, the product J Θ (x)•L v x yields a function vector whose components are sums of products of two sets of polynomials: those of degree up to q -1 from the Jacobian, and those of degree one from the linear transformation. The highest possible degree of any term in the product is q, as it is the sum of a polynomial of degree q -1 and a polynomial of degree one.Since Θ(x) includes all polynomials up to degree p, it includes all the terms generated by the productJ Θ (x) • L v x.Therefore, the components of J Θ (x)L v x ∈ R p can be written as linear combinations of the terms in Θ(x).C.1 Equivariant SINDy for Linear Symmetries C.1.1 Supplementary Results for Section 5.1In addition to Table 1 in the main text, Table 3 shows 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 error 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 a 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 L 2 error at each individual data point with high measurement noise. Figure 4 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 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.Our method can also be applied to higher-dimensional dynamical systems with symmetries. To demonstrate this, we consider the SEIR equations Abou-Ismail (2020) which model the evolution of pandemics by the number of susceptible (S), exposed (E), infected (I), and recovered (R) individuals.It is a 4-dimensional ODE system with quadratic terms. We use the following equations as the target of equation discovery:EQUATIONCompared 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 a moderate 5% noise to all trajectories.We assume an intuitive symmetry of scaling the number of recovered individuals proportional to the total population: v = (S + E + I + R)∂ R . We solve the constraint with regard to this symmetry as in Section 4.1. Figure 5 shows the reduced parameter space (from 60D to 34D). Table 4 shows the discovery results of SINDy and our EquivSINDy-c in terms of success probabilities and parameter estimation error 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. 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x) 1 x1 x2 x3 x4 x2 1 x1x2x1x3x1x4 x2 2 x2x3x2x4 x2 3 x3x4 x2 4 h1(x) h2(x) h3(x) h4(x)In addition to Table 2 in the main text, Figure 6 and Figure 7 show 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 to infinity. For the Lotka-Volterra system in Figure 6 , we observe that the prediction error grows slowest when using the discovered equations from our methods based on 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 a 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 a fluctuating error curve (after averaging). A potential reason is that it may discover some equations that are far from accurate but still lead to a periodic trajectory, similar to the true trajectory of the Sel'Kov glycolytic oscillator. The trajectory may have a large "eccentricity", so the error grows fast during a certain time and decreases rapidly afterward. Besides, we also visualize the simulated trajectories of some learned equations in Figure 8 . This provides a 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 8a , where predicted trajectories from most of the learned equations stay close to ground truth, but a few equations lead to diverging predictions due to the incorrect modeling of the exponential terms.Regularization Options. As mentioned in Section 4.2, we can introduce different forms of symmetry regularizations based on the four equations in Theorem 3.3 , listed as follows. The subscripts are short for Infinitesimal/Finite Group actions and Infinitesimal/Finite-time flow map of the ODE.EQUATIONEQUATIONEQUATIONEQUATIONwhere g = exp(ϵv) is computed for some fixed scale parameter ϵ. In other words, the summation is always performed over the learned infinitesimal generators (v ∈ B(g) and the group elements g in the above formulas are dependent on v.We also note that the denominator of (36) has a somewhat different form from others, i.e. the difference between f (gx) and f (x). If we use f (gx) or f (x) 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 ϵ → 0, τ → 0. The formulas inside the norms in the numerator and denominator in (35) both become O(ϵτ ), while f (gx) = f (exp(ϵv)x) = f (x + O(ϵ)) = x + O(ϵ) + O(τ ) + O(ϵτ ), whichdoes not match the numerator.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 disadvantages depending on the application, e.g. the representation of symmetry and the algorithm for equation discovery.In practice, we use only one of the losses (35, 36, 37, 38) 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 v therefore involves computing second-order derivatives of the network and amplifies the numerical error. On the other hand, if we have exact knowledge about the symmetry of the system in terms of its infinitesimal generators v, (38) becomes the most efficient way to incorporate symmetry loss, without the need of computing group elements based on the Lie algebra generators or integrating the learned ODE to its flow map. Besides, ( 37) and ( 38) have a 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) , it 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.the system towards a long horizon with these equations, the error will accumulate rapidly, as shown in Figure 3 . Also, the discovered equations from multiple runs are essentially different, making it difficult to select the best one to describe the dynamics.Compared to SINDyAE, symmetry-based 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 to satisfy the symmetry. The resulting equations can take slightly different forms, as shown in Table 7 . On the other hand, our method constantly discovers equations in the form ż1 = -az 2 , ż2 = bz 1 , where the numerical relationship between a and b depends on the discovered symmetries.We compare our method with the Time-Reversal 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 to encourage time-reversal symmetry. We evaluate TRS-ODEN as well as other baselines in Huh et al. (2020), i.e. Neural ODE (ODEN) and Hamiltonian Neural ODE (HODEN), on the damped oscillator system (13) and the Lotka-Volterra 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. Figure 9 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 our discovered equations with Lie symmetry regularization. Notably, in the Lotka-Volterra system, our method does not rely on any prior knowledge (e.g. time-reversal symmetry or 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.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 ∆t = 0.2. The initial conditions are sampled by sampling the distance to origin r ∼ U[0.5, 2] and the polar angle θ ∼ U[0, 2π] and then converting to (x 1 , x 2 ) = (r cos θ, r sin θ). The numerical integration is performed with a step size of 0.002. We apply white noise with noise level σ R = 20% to each dimension of the observation.Growth. We generate 100 trajectories for training, 20 for validation, and 20 for testing, with T = 100 timesteps sampled at a fixed rate ∆t = 0.02. The initial conditions are sampled uniformly from [0.2, 1] × [0.2, 1]. 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 xi = x i (1 + ϵ i ), where x i is the true state from numerical integration and ϵ i is sampled from a 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 a fixed rate ∆t = 0.002. For initial conditions, we first sample the population densities (p 1 , p 2 ) uniformly from [0, 1] × [0, 1]. Then they are converted to the logarithm population densities in the canonical form of the equation, i.e. x 1 = log(p 1 ) and x 2 = log(p 2 ). Also, we ensure that the Hamiltonian of the system H = exp(x 1 ) -x 1 + 1.333 exp(x 2 ) -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 a step size of 0.002. We apply white noise with noise level σ R = 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 ∆t = 0.002. The initial conditions are sampled uniformly from [0.5, 1] × [0.5, 1]. The numerical integration is performed with a step size of 0.002. We apply white noise with noise level σ R = 20% to each dimension of the observation. We note that our experimental results on this system are somewhat different from the results in Qian et al. (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 a result, the metrics such as success probabilities and parameter estimation errors appear worse than in Qian et al. (2022) . However, we still ensure fair comparison as we use the same generated dataset with our configurations for all the methods involved in our experiment.Reaction-Diffusion System. The system is governed by the PDEu t =(1 -(u 2 + v 2 ))u + (u 2 + v 2 )v + 0.1(u xx + u yy ), v t = -(u 2 + v 2 )u + (1 -(u 2 + v 2 ))v + 0.1(u xx + u yy ).We use the code from SINDy Autoencoder4 (Champion et al., 2019) to generate the data. The 2D space of u and v is discretized into a 100 × 100 grid. The system is simulated from one initial value for T = 6000 timesteps with step size ∆t = 0.05. Then, we add random Gaussian noises with standard deviation 10 -6 to each pixel and at each timestep. We use the timesteps t ∈ [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 ∈ [4800, 5400) to calculate the error.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-Volterra system where we also include the exponential terms. We apply sequential thresholding (Brunton et al., 2016) to enforce parsimony in the equations. The threshold is set to 0.05 for the damped oscillator and the growth system, 0.075 for the glycolytic oscillator, and 0.15 for the Lotka-Volterra system. The same threshold is applied to all methods including SINDy, WSINDy and EquivSINDy.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 no longer 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.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 on D-CODE with their official codebase5 . This codebase is developed based on 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.• If the contribution is a 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 is 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 to 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 a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • 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 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Section 5 discusses 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 of the paper to a 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, or as supplemental material.38th Conference on Neural Information Processing Systems (NeurIPS 2024).Note that v(f (x)) means evaluating the vector field at f (x) ∈ X. It is different from v(f )(x), which is the Lie derivative of f at the point x.Strictly speaking, the reverse only holds when ∆ satisfy certain additional local solvability conditions.https://github.com/kpchamp/SindyAutoencoders/tree/master/rd_solverhttps://github.com/ZhaozhiQIAN/D-CODE-ICLR-2022 


Decision: 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.

Reviews: review: paper_summary: paper_summary: This paper focuses on incorporating symmetries to equation discovery pipelines for ODEs. 

If 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. 

When 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.

When the symmetry is not known, the paper uses LaLiGAN, a symmetry discovery framework, to learn the corresponding regularization.
strength_weakness: strength_weakness: 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.
 
In 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.
questions: questions: What are the limitations of the approach in the known non-linear symmetries case?

What are the limitation of the approach in the unknown symmetries case?
limitations: limitations: The authors could comment on the technical limitations of their work (e.g. dimensionality, ODE order, sample complexity, etc).
paper_summary: paper_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
strength_weakness: strength_weakness: 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.

[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.
questions: questions: 1. Why set 5% and 20% noise respectively in the experiment in Section 5.1? What happens when the error scale is smaller?

2. For equation discovery with unknown symmetry, if we assume that the symmetry is unknown in the first two examples, what will be the result? 

3. The paper only considers two-dimensional situations. What are the results for high-dimensional examples?
limitations: limitations: Yes
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.  

Some 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.
questions: questions: Is there any existing theory on how well learned the symmetries can be from the data?

How does the computational cost of the method compare? I'm particularly interested in the increase in cost associated with the vector Jacobian product computation.

Is there any intuition why the RMSE of W-SINDy is better than the proposed method in the case of successful identification?
limitations: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: **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.

[1] Yang et al. Latent Space Symmetry Discovery. arXiv 2023.

***

**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.

***


There are some typos, e.g., the denominator of (35), and the summation spaces ($v$ instead of $g$) of (35) and (36)
questions: 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.

[2] Czarnecki et al. Sobolev Training for Neural Networks. NeurIPS 2017.
limitations: 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.
paper_summary: paper_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.
For this they make use of Lie groups.
strength_weakness: strength_weakness: _Weaknesses_
- 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.
- 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.


_Minor complaints/typos_
- A small arrow/text indicating the start "Start" in the flowchart on the left side of Figure 1 would make orientation slightly easier.
- Section 4 was harder to read, compared to section 3
- 261: With RK4 you mean Runge-Kutta 4 methods?
- Figure 2 could be bigger to improve readability.



To 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.
The most conflicting point for me are the results of Table 2, which I hope to have clarified during the discussion.
questions: questions: I want to highlight a major question before all others.
With respect to the results in Table 2:
- 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?
- 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?

Other questions:

- 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?
- 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$?
- For Theorem 3.3 I have a few minor questions:
	- What is $\epsilon$  in the $g = \exp(\epsilon v)$ term? 
	- Why is $g$ defined this way? Would it not be possible to calculate the Jacobian otherwise?
	- Is it guaranteed that $g \in G$? Does that only hold because $G$ is a symmetry group? Or generally?
	- 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?
- What is $L_v$ in Proposition 4.1? Is it a linear transformation? 
- 164/165: I don't understand what $c$ is and what its purpose is
- 180: The set of polynomials doesn't include products $x_1 x_2$ like in the example for SINDY in line 163, correct?
- 183: Is $W$ vectorized with respect to rows $\{[w_{11} ... w_{1n}], ...\}$  or columns $\{[w_{11} ... w_{n1}], ...\}$?
- Where does the factorization in Equation (10) come from? Are there established methods that easily generate this? What do its components mean?
limitations: limitations: In my opinion the only limitation of this work lies in the experiments in 5.3, as I stated before.
I 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.

Comments: comment: Thank you for the response. I keep my positive score of the paper.
Thank you for the authors' thoughtful response. Most of my concerns have been resolved, so I will be raising the review score.
Thank you for your thoughtful response! We agree that the regularization for non-linear symmetry does not *explicitly* constrain the parameter space. However, despite its difference from how we handle linear cases, this regularization approach is still shown to be effective in several dynamical systems in Sec 5.3 and Table 4 from the rebuttal PDF. By encouraging the learned equation to conform to the symmetry, even if not perfectly, we are much more likely to discover the correct equation. We will leave the problem of building explicit constraints for non-linear symmetries to future works.
Many thanks to the authors for the response, I have read and thought about it. On the linear comment, most of interesting parts of the paper are devoted to the linear case, which makes sense since this can be more tractably analyzed. The nonlinear case and its relative error regularization loss appear to be an afterthought. This is mainly because this regularization does not provide a hard constraint that the linear case does; rather, the symmetry is only satisfied to the extent that the regularization term is minimized (similar to a PINNs type regularization loss). Whereas in the linear case, the function is explicitly constrained to the space where the symmetries are satisfied, which is a stronger inductive bias.
Thank you for your response. I believe the authors responded well to my questions, and this paper is a valuable contribution to the area of learning dynamical systems. I have improved my score.
Thank you for answering my questions and providing clarifications.
I happily updated my score.