Discovering Symbolic Differential Equations with Symmetry Invariants
Authors: Jianke Yang
Manu Bhat
Bryan Hu
Yadi Cao
Nima Dehmamy
Robin Walters
Rose Yu
Venue: Accepted by TMLR ()

Abstract: Discovering symbolic differential equations from data uncovers fundamental dynamical laws underlying complex systems. However, existing methods often struggle with the vast search space of equations and may produce equations that violate known physical laws. In this work, we address these problems by introducing the concept of \textit{symmetry invariants} in equation discovery. We leverage the fact that differential equations admitting a symmetry group can be expressed in terms of differential invariants of symmetry transformations. Thus, we propose to use these invariants as atomic entities in equation discovery, ensuring the discovered equations satisfy the specified symmetry. Our approach integrates seamlessly with existing equation discovery methods such as sparse regression and genetic programming, improving their accuracy and efficiency. We validate the proposed method through applications to various physical systems, such as Darcy flow and reaction-diffusion, demonstrating its ability to recover parsimonious and interpretable equations that respect the laws of physics.

Published in Transactions on Machine Learning Research (03/2026)
Discovering Symbolic Differential Equations with Symmetry Invariants
Jianke Yang
jiy065 edu
University of Colifornia, San Diego
Manu Bhat
mbhat edu
University of California, San Diego
Bryan Hu
brhu@uced. edu
University of California, San Diego
Yadi Cao
yac0660rcsdd edu
University of California, San Diego
Nima Dehmamy
nima.dehmamy@idm.com
IBM Research
Robin Walters
r. walterss@orttheeeeee
Northeastern University
Rose Yu
roseyuSucad. edu
University of California, San Diego
Reviewed OpenReview: /ttps://openrevieee net/ forum? id= 9t1dByYfPc
Abstract
Discovering symbolic differential equations from data uncovers fundamental dynamical laws underlying complex systems. However, existing methods often struggle with the vast scarch space of equations and may produce equations that violate known physical laws. In this work, we address these problems by introducing the concept of symmetry invariants in equation discovery. We leverage the fact that differential equations admitting symmetry group can be expressed terms of differential invariants of symmetry transformations. Thus, we propose using these invariants as alomic entities in equation discovery, ensuring the discovered equations satisfy the specified symmetry. Our approach integrates seamlessly with existing equation discovery methods such sparse regression and genetic programming, improving their accuracy and elliciency. We validate the proposed method through applications to various physical systems, such as Darcy flow and -nfction-diffusion, demonstrating its ability to recover parsimonious and interpretable equations that respect the laws of physics.
1 Introduction
Differential equations describe relationships between functions representing physical quantities and their derivatives. They are crucial in modeling wide range of phenomena, from fluid dynamics and electromagnctic fields to chemical reactions and biological processes, as they succinctly caplure the underlying principles governing the behavior of complex systems. The discovery of governing equations in symbolic forms from obscrvational data bridges the gap between raw data and fundamental understanding of physical systems. Unlike black-box machine learning models, symbolic equations provide interpretable insights into the structure and dynamics of the systems of interest. In this paper, aim to discover symbolic partial differential

---

Published in Transactions on Machine Learning Research (03/2026)
Variables 1, a. 3r Mazi Gly
Symmetry Invariants 1, Au
Symmetry S0(2)
de 0.1
Enforced Symmetry
r:00
[ '0 't]
0:022 [11 44
Base
Algorithm
Mgy
Sparse Regression
Genetic Programming
Symbolic Transformers
Our framework enforces symmetry in equation discovery by using symmetry invariants. We highlight three discovery algorithms in their original form (bottom row) and when constrained to only use symmetry invariants (top row). The colored circles visualize the predicted functions on circular domain and demonstrate that using symmetry invariants guarantoes symmetric output.
Figure 1:
equations (PDFs) in the form
F(x, u (n)) = 0,
(1)
where x denotes the independent variables, w(n) consists of the dependent variable u and all of its up-to-nth order partial derivatives.
While it has long been exclusive task for human experts to identify governing equations, symbolic regression (SR) has emerged as an increasingly popular approach to automate the discovery[] SR constructs expressions from a predefined set of atomic entities, such as variables, constants, and mathematical operators, and fils the expressions to data by numerical optimization. Common methods include sparse regression [Brnnton et al.
2016 Champion et al., 2019), genetic programming (Cranmer et al. 2019 2020
Cranmer 2023
neural
networks [Kamienny et al. )022), etc.
However, symbolic regrossion algorithms may fail due to the vastness of the search space or produce more complex, less interpretable equations that overfit the data. A widely adopted remedy to these challenges is to incorporate inductive biases derived from physical laws, such as symmetry and conserved quantities, into equation discovery algorithms. Implementing these physical constraints narrows the space for equations and expodites the scarch process, and it also rules out physically invalid or unnecessarily complex equations.
Among the various physical constraints, symmetry plays fundamental role in physical systems, governing their invariances under transformations such as rotations, translations, and scaling. Previous research has shown the benefit of incorporating symmetry in equation discovery, such reducing the dimensionality of the search space and promoting parsimony in the discovered equation (Yang et al. 2024). However, the scopes of existing works exploiting symmetry are limited in terms of the types of equations they can handle, the compatible base algorithms, etc. For example, Udrescu & Tegmark (2020) dcals with algebraic equations; Otto et al. 2023) deals with ODE systems; Yang et al. (2024) applies to sparse regression but not other SR algorithms.
In this paper, we propose a general procedure based on symmelry invariands to enforce the inductive bias of symmetry with minimal restrictions in the types of equations and SR algorithms. Specifically, we leverage the fact that differential equation can be written in terms of the invariants of symmetry transformations if it admits a certain symmetry group. Thus, instead of operating on the original variables, our method uses the symmetry invariants as the atomic entities symbolic regression, as depicted in Figure [ These
While some literature symbolic regression specifically for GP-based methods, we the interchangeably with equation discovery to refer all sigorithms for learning symholic equations.
2

---

Published in Transactions on Machine Learning Research (03/2026)
invariants encapsulate the essential information while automatically satisfying the symmetry constraints. Consequently, the discovered equations are guaranteed to preserve the specified symmetry. In summary, our main contributions are listed as follows:
We propose a general framework to enforce symmetry in differential equation discovery based on the theory of differential invariants.
Our approach can be easily integrated with existing symbolic regression methods, such as sparse regression and genetic programming, and improves their accuracy and efficiency for differential equation discovery.
We show that our symmetry-basee approach is robust in challenging setups in equation discovery, such as noisy data and imperfect symmetry.
Notations. Throughout the paper, subscripts are usually reserved for partial derivatives, e.g. U == Ou/Ot, and Maz = Pu/Oxs. Superscripts are used for indexing vector components or list elements. We use Einstein notation, where repeated indices are summed over. Matrices, vectors and scalars are denoted by capital, bold and regular letters, respyctively, e-g. W,w, w. These conventions may admit exceptions for clarity or context. See Table 2 for a full description of notations.
2 Related Works
Symbolic Regression. Given the datasct {(at C X Y, symbolic regression (SR) aims to model the function y = fx) by symbolic equation. A popular method for symbolic regression is genetic programming (GP) (Schmidt & Lipson 2009; Gaucel et al. 2014), which leverages evolutionary algorithms to explore the space of possible equations and has demonstrated success in uncovering governing laws in various scientific domains such as material science (Wang et al.
2019), climate modcling (Grundner et al. 2023), cosmology
CCanmee et al. 2020), etc. Various soflware have been developed for GP-based symbolic regression, e.g.
Eurega (Dubčáková, 2011) and PySR (Cranmer
2023
Another class of methods is sparse regression [Brunton et 2016), which assumes the function to be discovered can be written as a linear combination of predefined candidate functions and solves for the cocfficient matrix. It has also been extended to discover more general equations, such equations in latent variables (Champion et al. 2019) and PDEs (Rudy et al. 2017).
Neural networks have also shown their potential in symbolic regression. Martius & Lampert (016)); Sahoo ct al. (2018) represents a few carliest attempts, where they replace the activation functions in fully connected networks with math operators and functions, so the network itself translates to a symbolic formula. Other works represent mathematical expressions as sequences of tokens and train neural networks to predict the sequence given a dataset of input-output pairs. For example, Petersen al. (2019) trains an RNN with policy gradients to minimize the regression error. [Biggio et al. [2021). Kamienny et al. (2022) and Holt et al. (223)) pre-train an encoder-decodee network over a large amount of procedurally generated equations and query the pretrained model on a new dataset of input-output pairs at test time.
The aforementioned symbolic regression methods can be improved by incorporating specific domain knowledge. For example, AT Feynman (Udrescu & Tegmark 2020 Udrescu et al. 2020) uses properties like separability and compositionality to simplify the data. Cranmer et al. (2020) specifics the overall skeleton of the equation and fits each part with genetic programming independently. The goal of this paper falls into this category to use the knowledge of symmetry to reduce the search space of symbolic regression and improve its accuracy and efficiency.
Recently, Large Language Models (LLMs) have emerged as an alternative for SR, using pre-trained scientific priors to propose sequential hypothesis (Merler et al. 2024) or to guide genetic programming (Shojaee et al,, 2024), balancing the efficiency of domain knowledge with the robustness of evolutionary scarch. However, current LLM-based methods often rely on memorizing known equations rather than facilitating genuine discovery, and their guidance lacks interpretability, specifically, the reasoning behind their suggestions,
3

---

Published in Transactions on Machine Learning Research (03/2026)
evidenced by recent benchmark specially designed for LLM-SR (Shojace al. 2025 A recent effort sought. to improve interpretability by binding symbolic evolution with natural language explanations (Grayeli et al. 2024). However, this method relies on frontier LLMs to conduct the evolution of the natural language components, rendering the process itself opaque. These limitations highlight the need for approaches that enhance the controllability and explainability of the prior knowledge injected, ensuring more transparent and trustworthy discovery.
Discovering Differential Equations. While it remains in the scope of symbolic regression, the discovery of differential equations poses additional challenges because the derivatives are not diroctly observed from data. Building upon the aforementioned SINDy sparse regression (Brunton et al., 2016), Messenger & Bortz (2021a b) formulates an alternative optimization problem based on the variational form of differential equations and bypasses the need for derivative estimation. similar variational approach is also applied to genetic programming (Qian et al. 2022). Various other improvements have been made, including refined training procedure (Rao et al.) 2022), relaxed assumptions about the form of the equation (Kaheman et al.
2020), and the incorporation of physical priors (Xie et al., 2022 Bakarji al. 2022 Lee al. 2022 Messenger
ct al. 2024
PDE Learning and Surrogate Modeling. Beyond symbolic regression, there is broad line of work learning PDEs and their solution operators directly from data using non-symbolic surrogates. Neural operator methods aim to approximate noulinear operators mapping initial or boundary data to PDE solutions, providing highly expressive black-box solvers for families of PDEs (Lu et al.) 2021 Li et al., 2021). Gaussian processes and kernel methods have also been developed to solve nonlinear PDEs and related inverse problems (Chen ct al. 2021). A complementary body of work focuses on learning effective coarse- gaarse-griindd dynamics or closures from fine-scale simulations and spatiotemporal data, including reservior-computing and recurrent architectures Vlachas et al. 2020), machine-learned coarse-scale PDES from microscopic data (Lee et al. 2020), numerical bifurcation analysis and intrinsic-coordinat models inferred from simulators (Galaris et al,
2022 Floryan & Graham 2022), and multiscale frameworks that learn reduced stochastic or PDE models for
complex systems (Vlachas et al. 2022: Lee et al., 2023 Dietrich et al., 2023 Fabiani et al. 2024). These approaches typically prioritize predictive accuracy and efficient surrogate modeling over known symnollc structure. In contrast, our framework, along with other SR methods, seeks to recover interpretable closed-form PDEs whose terms can be inspected and analyzed, while still leveraging data-driven tools.
PDE Symmetry in Machine Learning. Symmetry is an important inductive bias in machine learning. In the context of learning differential equation systems, many works encourage symmetry in their models through data augmentation (Brandstetter et al. 2022), regularization terms Likhound-Sadege et al. 2023
Zhang et al. 2023 Dalton et al. 2024), and sell-supervised learning (Mialon et al. 2023). Strictly enforcing symmetry is also possible, but is often restricted to specific symmetries and systems (Wang et al., 2021 Gurevich et al. 2024). For more general symmctrics and physical systems, enforcing symmetry often requires additional assunptions on the form of equations, such as the linear combination form in sparse regression Otto et al., 2023 Yang al. 2024h. To the best of our knowledge, our work is the first attempt to strictly enforce general symmetries of differential equations for general symbolic regression methods. A more detailed discussion of the connections and differences between our work and other symmetry- based equation discovery methods is provided in A Appendix] D
3
Background
3.1
PDE Symmetry
This section introduces the basic concepts of partial differential equations and their symmetry. For a more thorough understanding of Lie point symmetry of PDEs, we refer the readers to Olver (1993).
Partial Diffcrential Equations. We consider PDEs in the form F(x, mn))) = 0, as given We restrict ourselves to single equation and a single dependent variable here, though generalization multiple equations and dependent variables is possible. We use x E X c RRP to denote all independent. variables. For example,
4

---

Published in Transactions on Machine Learning Research (03/2026)
x= (t, for system evolving in 1D space. Note that the bold refers to the collection of all independent variables while the regular x denotes the spatial variable. Then, u(x) EU C IR is the dependent variable; u(x) = (u, uay... denotes all up to nth-order partial derivatives of u; (x, u(")) E M(") C X x Utm), where M(n) is the nth order jet space of the total space X x U. M(") and u(") are also known as the ath-order prolongation of X x U and u, respectively.
Symmetry of a PDE. A point symmetry g is a local diffcomorphism on the total space = X x U:
g (x, u) = (*(x, u), ū(x; u)),
where x and à are functions on E. The action of g on the function u(x) is induced from ( by applying it to the graph of u X + U. Specifically, denote the domain of as C X and its graph as I,, ((x, u(x)) xx E @l).
The group element g transforms the graph r, as F. = g -u = {(x, ū) = 9 (x, u) (x, u) E Tu).
Since g transforms both independent and dependent variables, I' does not necessarily correspond to the graph of any single-valued function. Nevertheless, by suitably shrinking the domain Ω can ensure that the transformations close to the identity transform Ta to the graph of another function. This function with the transformed graph F. is then defined to be the transformed function of the original solution u, i.e. g u = à s.t. r. = F'.. The symmetry of the PDE I is then defined:
Definition 3.1. A symmetry group of F(x, u!")) = 0 is local group of transformations G acting on an open subset of the total space x U such that, for any solution u to F = and any g E G, the function = (g (x) is also a solution of = wherever
defined.
Infinitesimal Generators. Often, the symmetry group of PDE is continuous Lie group. In practice, needs to compute with infinitesimal generators of continuous symmetries, ie., vector fields. In more detail, we will write vector fields v E + TE on E = x U as
a
v= (x, u) dxi + o(x, u)
(3)
Ou
Any such vector field generates a one-parameter group of symmetries of the total space [exp(ev) : E IRJ. The symmetries arising from the exponentiation of a vector field moves a point in the lotal space along the directions given by the vector field. We will specify symmetries by vector fields in the following sections. For instance, v= vr0 yd, represents the rotation in (z, y)-pane; v= B corresponds to time translation.
To analyze the symmetry of PDEs, we must know how it transforms not only the variables, but also their derivatives accordingly. The group transformations derivatives are formalized by prolonged group actions and infinitesimal actions on the nth-order jet space, denoted g(*) and v(*), respectively. More details on prolonged group actions are discussed in Appendix A.2, with Figure 5 visualizing a simple example. To
introduce our method, it suffices to note that the prolongation of the vector field 3 can be described explicitly by ş and & and their derivatives via the prolongation formula
3.2 Symbolic Regression Algorithms
Given the data y C x, the objective of symbolic regression (SR) is to find a symbolic expression for the function y f(). Although this original formulation is for algebraic equations, it can be generalized
to differential equations like [ To discover a PDE from the dataset of its observed solutions on a grid Q, i.c. [(x u(x)) : x E [], we estimate the partial derivative terms and add them to the dataset: {(x, 2())) x E Ω). One of the variables in the variable set (x, u(") is used as the LHS of the equation, i.e. the role of the label y in symbolic regression, while other variables serve as features. The precise set of derivatives added to symbolic regression and the choice of the equation LHS requires prior knowledge speculations about the underlying system.
We briefly review two classes of SR algorithms: sparse regression (SINDy) and genetic programming (GP).
Sparse regression (Brunton et al. 2016) is specifically designed for discovering differential equations. It assumes the L.HS & of the equation is fixed term, e.g. f = ue, and the RHS of the equation can be written
5

---

Published in Transactions on Machine Learning Research (03/2026)
linear combination of m predefined functions @ with trainable coefficients R", ie.,
f(x, u (m))) = w 0 (xx ) 03 M(n) R.
(4)
The equation is found by solving for w that minimizes the objective I|- R + A/|w/lo: where L and R are obtained by evaluating and wisi on all data points and concatenating them into column vectors, and ||wloo regularizes the number of nonzero lerms. This formulation can be easily extended to q equations and dependent variables (q 1): ((x, u())) w (x, u())., W R
One problem with sparse regression is its restrictive assumptions about the form of equations. Many equations cannot be expressed in the form of (4) e.g. y - where a could be any constant. Also, the success of sparse regression relies on the proper choice of the function library {0i]. If any term in the true equation were not included, sparse regression would fail to identify the correct equation.
Genetic programming (GP) offers alternative solution for SR (Cranmer 2023 which can learn equations in more general forms. It represents each expression as a tree and instantiates a population of individual expressions. At each iteration, it samples a subset of expressions and selects one of them that best fits the data; the selected expression is then mutated by a random mutation, a crossover with another expression, or a constant optimization; the mutated expression replaces an expression in the population that does not fit the data well. The algorithm repeats this process to search for different combinations of variables, constants, and operators, and finally returns the "fittest" expression. GP can be less efficient than STNDy when the equation can be expressed in the form 4 due to its larger search space. However, we will show that it is a promising allernative to discover PDEs of generic forms, and our approach further boosts ils efficiency.
4 Symbolic Regression with Symmetry Invariants
Symmetry offers a natural inductive bias for the search space of symbolic regression in differential equations. It reduces the dimensionality of the space and encourages parsimony of the resulting equations. To enforce symmetry in PDE discovery, we aim to find the maximal set of equations admitting a given symmetry and search in that set with symbolic regression (SR) methods.
4.1 Differential Invariants and Symmetry Conditions
To achieve this, our general strategy is to replace the original variable set with a complete set of invariant functions of the given symmetry group. Since we consider PDEs containing partial derivatives, the invariant functions refer to the differential invariants defined as follows.
Definition 4.1 (Def. 2.51, Olver [1993). Let G be local group of transformations acting on X x U. Any E G gives a prolonged group action (r(")g g on the jet space (()) C X x U()), An nth order differential invariant of G is a smooth function 1 : (n) - R, such that for all E G and all (x, (()) E M(0) r(g(")) (x, m()))) = m)(x, (n)) whenever g(m). (x, (")) is defined.
In other words, differential invariants are functions of all variables and partial derivatives that remain invariant under prolonged group actions. Equivalently, if G generated by set of infinitesimal generators B = {va), then a function is a differential invariant of G iff (()) (n) = 0 for all v. E B. The following thcorem guarantees that any differential equation admitting a symmetry group can be expressed solely in terms of the group invariants.
Theorem 4.2 (Prop. 2.56, Olver 1993)). Let G be local group of transformations acting on X x U. Let (x, (u.,)),..,?"/(x } be a complete set of functionally independent ntthorddr differential invariants of G. An nth-order differential equation {1] admits G as symmetry group if and only if it is equivalent to equation of the form F(n) .?..,?"** O.
Consequently, SR with a complete set of invariants precisely searches within the space of all symmetric differential equations and automatically excludes equations violating the specified symmetry.
Our strategy of using differential invariants applics broadly to various equation discovery algorithms. For instance, in sparse regression, construct the function library using invariants rather than variables
6

---

Published in Transactions on Machine Learning Research (03/2026)
and derivatives. Similarly, in genetic programming, the variable set can be redefined to include only invariant functions. In each case, the key benefit. is the same: the search space is restricted to symmetry-respecting equations by construction. The reduced complexity of the equation search also leads to increased accuracy and efficiency.
Next, we describe how to construct a complete set of differential invariants (Section 4.2 and how incorporate them into specific SR algorithms (Section 4.3
4.2 Constructing Complete Set of Invariants
Despite the simplicity of our strategy, we still need concrete method for computing the invariants. In this subsection, we provide general guideline to construct complete set of differential invariants up to required order given the group action.
By definition of differential invariants, we look for functions n7(x, u(")) satisfying (()) = 0 given a prolonged vector field v(m), This is first-order linear PDE that can be solved by the method of characteristics. However, in practice, E = X x U ~ R' x R, there are (pta-1) partial derivatives of the independent variable u of order exactly n. Therefore, as n grows, quickly becomes impractical to solve directly for nth-order differential invariants. The higher-order differential invariants, if necessary, can be computed recursively from lower-order ones by the following result:
Let ,/,,,"" ,""..... .pp be any differential invoriants of G whose horizontal Jacobian J= [Dr'] is non non-degenerate on an open subset Ω C M()). If there are a mazimal number of independent, strictly nth-order differential invoriants C).... S On = -(pt (p+n-1), then the following set contains a complete set of independent, striclly (n+ 1o+1)th-ordee differential invariants defined on S2:
Proposition 4.3. Let G be local group of tvansformations acting on
det (ee tk E E [ggn],
(5)
where j,j E [p] are matria indices, D; denoles the lotal derivalive w.r.l i-th independent variable and 1 (/ = |n "...,,"-,, + ..... .").
In practice, first solve for pr v(n) 0 to obtain sufficient number of lower-order invariants required in Proposition
4.3 and then construct complete sets of invariants of arbitrary orders. Notably, while in theory our framework operates on any complete set of differential invariants, the invariants computed this way may be algebraically complicated and poorly scaled, leading to difficultics in SR optimization. In practice, we start from such a complete set of differential invariants and then deliberately convert them into simpler, physically interpretable invariant functions (such as Laplacians for rotational symmetry) as the feature set for SR. Then, we evaluate invariants on the dataset only where they are well-defined. If necessary, we shrink the domain and filter out data points that cause singularity (e.g., where the denominator of an invariant function vanishes). In Appendix A.4, we provide two examples of different symmetry groups and their differential invariants. Those results will also be used in our experiments.
4.3 Implementation in SR Algorithms
Our symmetry principle characterizes a subspace of all equations with given symmetry. Generally, this subspace partially overlaps with the hypothesis spaces of SR algorithms, conceptually visualized in Figure 2 As in Theorem 4.2 PDEs with symmetry can be expressed as implicit functions of all differential invariants. However, symbolic regression methods typically learn erplicit functions mapping features to labels. Some algorithms, such SINDy, impose even stronger constraints equation forms. Therefore, adaptation is needed to implement our strategy of using differential invariants in specific symbolic regression algorithms. In general, it requires more effort to adapt our method to base algorithms with more structured and nontrivial hypotheses about possible equation forms. Below, we discuss in detail how adapt our method to two common classes of base SR algorithms.
All equations: alx.
SR: Flec.)
SINDY: werz
Equations w symmelty gi )
Symmetry SINDY
Symmetry n SR
Figure 2: Venn diagram of hypothe- sis spaces from base SR methods and our symmetry principle.
7

---

Published in Transactions on Machine Learning Research (03/2026)
General explicit SR We start with general SR methods that learn erplicit function y= f(z) without additional assumptions about the form of J, e.g.; genetic programming and symbolic transformer. When learning the equation with differential invariants, we do not know which one of them should be used as the LHS of the equation, i.c. the label in symbolic regression. Thus, we fit an equation for cach invariant as LHS and choose the equation with the lowest data error, as described in Algorithm | to select the best equation since the scales of LHS terms differ.
We use relative error
Algorithm General explicit SR for differential equations with symmetry invariants
Require: PDE order n, datasct {z' = (x', (u(m)) M( No /ill base SR algorithm S: (X,)) + y f(x),
infinitesimal generators of the symmetry group B = [va).
Ensure: A PDE admitting the given symmetry group. Compute the symmetry invariants of B up to nth-order: 17", ..... ,77 {Prop. 4.3 Evaluate the invariant functions on the dataset: ry = (z'), for E [K], i [Np).
Initialize a list of candidate equations and their risks: E= D. for & in 1 K do
Use the kth invariant as label and the rest as features: y m X 17- k,;
Run S(X. y) and get candidate equation n/ = f*(n k) Evaluate L =
- f (X)l///llll and set E[k] = (f*, (*).
end for
Choose the equation in with the lowest error: h = arg min, B[j][2].
return 17" = (n *). {Optionally, expand all np² in terms of original variables z.}
Sparse regression SINDy assumes linear equation form 4 Generally, its function library differs from the set of differential invariants. Also, SINDy fixes LHS term, while we do not single out invariant the LHS of the equation when constructing the set of invariants.
Assume we are provided the SINDY configuration, iee the LHS term l and the function library {0'). To implement sparse regression with symmelry invariants, we assign an invariant 1/* that symbolically depends on f, Le. ônz / * O, as the LHS for the equation in terms of symmetry invariants. The remaining invariants are included on the RHS, where they serve as inputs of the original SINDy library functions. In other words, the equation form is 1/* w. (0) (n) Similar to Algorithm we can expand all v variables to obtain the equation in original jet variables.
The above approach optimizes an unconstrained coefficient vector w for functions of symmetry invariants. Alternatively, we can use the original SINDy equation form 4 and implement the symmetry constraint as a constraint on the coefficient w, as demonstrated in the following theorem. Here, we generalize the setup to multiple dependent variables and equations.
Proposition 4.4. Let e(x, ((n))= O u()) be system of q differential equations admitting a symmetry group G, where x E RR, u E R2, 0 E IRm Assume there erist some nilh-order invariants of G, 1lo 1:g and 17 i,t.K, s.t. (1) the system of equations can be expressed as тo = w'''(m),, where no = [no 1:91 and n = =b and (2) Th = gi and (#')²= SS for some functions ('()) and constant tensors W T and S. Then, the space of all possible W is lineur subspace of Raxm.
Intuitively, the conditions above state that the equations can be expressed as a linear combination of invariant. terms, similar to the form 4 w.r.t original jet variables. Also, every invariant term in no and 0' ('(1)) is already encoded in the original library O. In practice, we need to choose a suitable set of invariants according to the SINDy configuration to meet these conditions. For example, it is a common SINDy setup where O contains all monomials on Mn)) up to degree d. In this case, any set of invariants where each invariant is a polynomial on (x) up to degree d satisfies these conditions.
The proof of Proposition| 4.4 is provided in Appendix AA.5] where we explicitly identify the basis of the linear subspace for W. Then, we can use this basis to build a SINDy model for PDEs with the corresponding symmetry, similar to how EMLP (Finzi et al.) (021b) constructs equivariant linear layers and how Yang et al. (2024) construets equivariant SINDy for ODFs. Specifically, if the constrained subspace has basis
8

---

Published in Transactions on Machine Learning Research (03/2026)
Q R=xxx" where is the subspace dimension, write W = QukB'. We then fix Q and solve 8 R" using the same least square optimizer in SINDy, effectively reducing the parameter space dimension from X m r.
In practice, we observe that the basis Q obtained from the constructive proof of Proposition 4.4 is not sparse. The lack of sparsity can pose problem when we perform sequential thresholding in sparse regression. Specifically, in SINDy, the entries in W that are close zero are filtered out at the end of each iteration, which serves a proxy to the sparsity-prmmoting Ln regularization. Since we fix Q and only optimize B, straightforward modification to the sequential thresholding procedure is to threshold the entries in 3 instead of those in W. However, if Q is dense, even a sparse vector 3 can lead to dense W, which contradicts the purpose of sparse regression. To address this issue, we apply a Sparse PCA to Q to obtain a sparsified basis. More implementation details and examples of the computed basis Q can be found in Appendix B.2
One notable advantage of converting symmetry into lincar constraints with Proposition 4.4 that allows to keep track of the original SINDy parameters W during optimization. This enables straight.orward integration of symmetry constraints to variants of SINDy, e.g. Weak SINDy (Messenger & Bortz 2021a for noisy data. For Wik = =Qkk while directly optimize B, we can still easily compute the objective of Weak SINDy which explicitly depends on W. In comparison, if we use the raw invariant terms for regression, e.g. the equations lake the form Yo- = (')'(n), it is challenging to formulate the objective of Weak SINDy with respect to W', In Appendix B.3 we provide more details on how to implement this constrannt
approach on Weak SINDy.
4.4 Constraint Relaxation for Systems with Imperfect Symmetry
Our approach discovers PDEs assuming perfect symmetry. However, it is common in reality that system shows imperfect symmetry due to external forces, boundary conditions, etc. (Wang et al.) 2022). In such cases, the previous method cannot identify any grmmettyybreaking factors.
To address this, we propose to relax the symmetry constraints by allowing symmetry-breaking terms to appear in the equation, but at a higher "cost". We implement this idea in sparse regression, where the equation has a linear structure = We. We adopt the technique from Residual Pathway Prior (RPP) (Finzi et al. 2021a), which is originally developed for equivariant linear layers in neural networks. Specifically, let Q be the basis of the parameter subspace that preserves symmetry and P be the orthogonal complement of Q. Instead of parameterizing W in this subspace, we define W = A + B where A = Q 3 and Bik pijky', where B and y are learnable parameters, and place a stronger regularization on than on B. While the model still favors equations in the symmetry subspace spanned by Q, symmetry-breaking components in P can appear fits the data well.
5
Experiments
5.1 Datasets and Their Symmetries
We consider the following PDE systems, which cover different challenges in PDE discovery, such high- order derivatives, generic equation form, multiple dependent variables and equations, noisy dataset, and imperfect symmetry. The datasets are generated by simulating the ground truth equation from specified initial conditions, with detailed procedures described in Appendix E.1
Boussinesq Equation. Consider the Boussinesq equation describing the unidirectional propagation of solitary wave in shallow water (Newell 1985):
+ UUlz + u + Uazza
(6)
This equation has scaling symmetry V 2td, rr,-- 2u0,, and the translation symmetries in space and time. The differential invariants are given by Ma.9) = ( ( (22+0+28)/3 where a and / are the orders of partial derivatives in æ and t, respectively. To discover the 4th-order equation, we compute all Ma,,) for 0 + B & 4, except for 71,00) 1 which is a constant.
9

---

Published in Transactions on Machine Learning Research (03/2026)
Darcy Flow. The following PDE describes the steady state of 2D Darey flow (Takamoto et al. 2022 with
spatially varying viscosity a(z, 3) 4(x +v") and a constant force term f(x) =1
4(x vu)
(7)
This equation admits an S0(2) rotation symmetry V yd, xdyy A detailed calculation of the differential
invariants of this group can be found in Example A.5 In experiment, use the following complete set of
2nd-order invariants: (a +2²), u, M, yM 2 yily Mex + Myy: u + 2u uu 2 x² Uax + uyu 2nyuwy
Reaction-Diffusion. We consider the following system of PDEs from Champion et al. 2019
0 = d (1 (1- - v²)u (u² + v²)y 2 = d2 v - (u² + vJ)u + (1 u² - v²)u
(8)
In the default setup, we use d1 d2 0.1. The system then exhibits rotational symmetry in the phase space: ud, vθ,. The ordinary invariants (functions of variables, not derivatives) are (1,2,,, y, u² + v²). The higher-order invariants are {u u, u- UnJ. where vu,")² and p is any multi-index of t, æ and
We also consider the following cases where the rotation symmetry is broken due to different factors:
Unequal diffusivities We use different diffusion coefficients for the two components: dj = 0.1, d2 = 0.1 + This can happen, for example, when two chemical species described by the equation diffuse at different rates due to molecular size, charge, or solvent interactions.
External forcing The ground truth equation (8) is modified by adding -c to the RHS of u, and -eu to the RHS of er. This can reflect a weak parametric forcing on the system.
5.2 Methods and Evaluation Criteria
We consider three classes of algorithms for equation discovery: sparse regression (PySINDy, de Silva et al. (2020); Kaptanoglu et al. )), genetic programming (PySR, Cranmer (2023 ), and a pretrained symbolic transformer (E2E, Kamienny et al. (2022)). For each class, we compare the original algorithm using the regular jet space variables (i.e. (x, m()) and our method using symmetry invariants. Our method will be referenced SI (Symmetry Invariants) in the results.
To evaluate equation discovery algorithm, run it 100 times with randomly sampled data subsets and randomly initialized models if applicable. We record its success probubilily (SP) of discovering the correct equation. Specifically, we expand the ground truth equation into c(f = O, where of are nonzero coefficients, z denotes the variables involved in the algorithm, i.c., original jet variables (x, u()) for basclines and symmetry invariants for our method, and I' are functions of z. Also, the discovered equation is expanded as Σ, (z) 0, where at + The discovered equation considered correct if all the terms with nonzero coefficients match the ground truth, i.e., (/'} = (f*. We also report the prediction error (PE), which measures how well the discovered equation fits the data. For evolution equations with time derivatives on the LHS, we simulate each discovered equation from an initial condition and mcasure its difference from the ground truth solution a specific timestep in terms of root mean square error (RMSE). Otherwise, we just report the RMSE of the discovered equation evaluated all test data points.
10

---

Published in Transactions on Machine Learning Research (03/2026)
5.3 Results Clean Data with Perfect Symmetry
Table 1: Equation discovery results on clean data. C, standing for complexity, refers to the effective parameter space dimension in sparse regression and the number of variables in GP GP/Transformer SP and PE stands for success probability and prediction error, as explained in Section 5.2 The entries *-* suggest that the method does not apply to the specific PDE system, or the result is not meaningful. The arrows f mcan higher/lower metrics are better. Confidence intervals for PE are reported in Table
Boussinesq SP f
66
Darcy flow
7 PE
Reaction-diffusion
(8
Method
C
+
PE 4
C
+
SP
+
C
SP 1
PE +
Sparse
PySINDy
15
0.00
0.373
38
0.53
0.021
Regression
SI
15
1.00
0.098
28
0.54
0.008
Genetic
PySR
17
0.90
0.098
8
0.00
0.114
17
0.00
Programming
SI
14
1.00
0.098
7
0.79
0.051
16
0.81
0.023
E2E
10
0.53
0.132
8
0.00
17
0.00
Transformer
SI
0.85
0.104
7
0.00
16
0.00
Table 1 summarizes the performance of all methods the three PDE systems. For prediction errors (PE), we report the median, instead of the mean + standard deviation, of 100 runs for each algorithm. This is because some incorrectly discovered equations yield large prediction errors, making the mean and the standard deviation less meaningful. We additionally report the confidence intervals for prediction errors by (25% quantile, 75% quantile] in Table 3 Comparisons are made within each class of methods. Generally, using symmetry invariants reduces the complexity, defined as the effective parameter space dimension sparse regression and the number of variables in GP /rrnnsfommrr, of equation discovery, and improves the chance of finding the correct equations compared to the baselines.
Specifically, in sparse regression, our method (SI) using symmetry invariants is only slightly better than PySINDy in the reaction-diflusion system, but constantly succeeds in the Boussinesq equation where PySINDy fails. The failure PySINDy is because the ut term in (6) is not supported by its function library, showing that SINDy's success relies heavily on the choice of function library. This exclusion of function terms also explains why the reported complexity (C) of PySINDy the same as that our mothod (SI), which would otherwise be smaller because symmetry reduces the hypothesis space of equations. On the other hand, by enforcing the equation to be expressed in invariants, our method automatically identifies the proper function library. Appendix C.2 provides results for other variants of sparse regression, where we modify the implementation of PySINDy to include different terms in its library.
For GP-based methods, Table displays the results with a fixed number of GP iterations for cach datasct. We also experiment with different numbers of iterations in Figure 3 showing that our method can identify the correct equations with fewer iterations and is considered more efficient.
Boussinesq GP Success Probabiitty
Darcy GP Success Probability
R-D GP Success Probability
SI (ours)
(ours)
Py5R
PySR
Prabatiity
2.5
0.2
500005 Probabiiityy
2.6
2.4
0.4
2.4
Si (ours) PySR
8.9
0.0
8.0
35
15
200
203
400
no
Iterattons
Figure 3: Success Probabilities of GP-based methods on different systems. Our method with symmetry invariants can discover the correct equations with fewer iterations.
11

---

Published in Transactions on Machine Learning Research (03/2026)
On the other hand, the pretrained symbolic transformer fails two of the three datasets. We conjecture this is because the data distribution from PDE solutions greatly differs from its pretraining dataset. However, the symbolic transformer can discover the Boussinesq equation correctly, where using symmetry invariants leads to a much higher success probability.
5.4 Results on Noisy Data and Imperfect Symmetry
We test the robustness of our method under two challenging scenarios: (1) noise in observed data, and (2) PDE with imperfect symmetry.
R-D w/ noisy data
R-D unequal diffusivities
R-D w/ external forcing
0.3
SI (autotune)
SI (defa ult)
0.5
WSINDY fautoture)
Success Probability
WSINDy (default)
Success probability
Success probability
0.0
0.4
0.2
6.4
St-relaxed
6.1
Si-re axed
SI
WSINDD
WSINDY
Perfect symmetyy
Perfect symmetry
0.0
02-
2.5
0.02
0.02
-.21
5
0.02
0.03
0.1
tis
Noise Level (96)
Symmetry breaking parameter
Symmetry breaking parameter
Figure 4: Success probabilities of sparse regression methods on the reaction-diffusion system with noisy data (left), unequal diffusivities (center) and external forcing (right). Under noisy data, our method (SI) consistently outperforms WSINDy under the same test. function configuration. For systems with imperfect. symmetry, strictly enforcing symmetry (SI) hurts performance, but a relaxed symmetry constraint (SI-relaxed, introduced in Section 4.4h still better than no inductive bias (WSINDy).
In the first experiment, add different levels of white noise to the simulated solution of the reaction-diffusion system. Since the derivatives estimated by finite difference are inaccurate with the noisy solution, we use WSINDy (Messenger& & Bortz 2021a), which does not require derivative estimation. The success probabilities of our method (SI) and WSINDy are shown in Figure 4 (left). We adopt two different strategies for choosing the test function parameters in WSINDy, including the number of test functions, half-lengths of the square subdomains in each spatiotemporal direction, and the polynomial degree of test functions: (1) following Messenger & Bortz (2021a), the parameters aulo-luned based on the noise characteristic of the data; and (2) the default parameters in PySINDy used. In both WSINDy setups,
method consistently achieves higher success probability at different noise levels. Notably, when the noise level is high, our symmetry-conssrrinee model performs better with the default parameter setup, which, as we observe, uses fewer test. functions and lower polynomial degrees than the auto-tuned setup. We comment that choosing test functions and related hyperparameters is known to be a challenging problem (Bortz et al. 2023 Tran & Bortz 2025), and we leave further investigation of this phenomenon to future work.
In the second experiment,
simulate the two variants of { (unequal diffusivities and external forcing) with different values for the symmetry-breaking parameter and add 2% noise to the numerical solutions. We compare three models: (1) our model with strictly enforced symmetry (ST), (2) our model with relaxed symmetry (SI-relaxed) introduced in Section 4.4 and (3) WSINDy with default PySINDy parameters as the bascline. The results for the two systoms with symmetry breaking are shown in Figure ! (center & right). As expected, ST has a much lower success probability when the symmetryy-reakee Mcanwhile, SI-relaxed remains highly competitive. It also has clear ad advantage over baseline SINDy, showing that even if the inductive bias of symmetry is slightly inaccurate, our model with relaxed constraints is still better than model without any knowledge of symmetry. We also include the results the same systems with auto-tuned WSINDy parameters in Appendix C.3
factor becomes significant.
More comprehensive results, e.g. variant sparse regression models, comparison with an additional D-CIPHER Kacpızyk et al. 2023) baseline, equations with higher spatial dimensions and larger symmetry groups, discovered equation samples, are provided in Appendix
12

---

Published in Transactions on Machine Learning Research (03/2026)
6 Discussion
We propose to enforce symmetry in symbolic regression algorithms for discovering PDEs by using differential invariants of the symmetry group as the variable set. We implement this general strategy in different classes of algorithms and observe improved accuracy, efficiency and robust.nss of equation discovery, especially in challenging scenarios such as noisy data and imperfect symmetry.
It should be noted that our method assumes the symmetry group is already given. This assumption aligns with common practice: physicists often begin by hypothesizing the symmetries of a system and seek governing equations allowed by those symmetries. However, unknown, and will produce incorrect results with misspecified symmetry. This can be potentially addressed by incorporating automated symmetry discovery methods for differential equations Yang et al. 2024 Ko
current framework camnot be applied if symmetry is
et al. 2024); which we leave for future work.
Another caveat our method is the calculation of differential invariants. While solving for (n) = 0 and applying the formula (5) is easy with any symbolic computation package, the resulting differential invariants may be complicated and require ad-hoc adjustment for better interpretability and compatibility with specific algorithm implementations (e.g. conditions in Proposition 4.4p. Fortunately, this only requires one-time effort. Once we have derived the invariants for symmetry group, the results can be reused for any equation admitting the same symmetry.
Broader Impact Statement
The method in this paper can potentially be used to expedite the process of discovering governing equations from data and aid researchers in other scientific domains. Equally important, equations inferred from imperfect or biased data may appear authoritative yet embed systematic errors. Thorough validation checks, uncertainty quantification, and domain-expert review protocols for the discovered equations are essential.
References
Tara Akhound-Sadegh, Laurence Perreault-Levasseun Johannes Brandstetter, Max Welling, and Siamak Ravanbakhsh. Lie point symmetry and physics-informed networks. Advances in Neural Information Processing Systems, 36:42468 42481, 2023.
Joseph Bakarji, Jared Callaham, Steven L. Brunton, and Nathan Kutz. Dimensionally consistent learning with buckingham pi. Nature Computational Science, 2:834-844 2022. ISSN 2662-8457. doi: 10.1038/
4445585-022-033555
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In Marina Mcila and Tong Zhang (eds.). Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 936-945.
PMLR, 18 24 Jul 2021.
David M Bortz, Daniel A Messenger, and Vanja Dukic Direct estimation of parameters in ode models using wendy: Weak-form estimation of nonlinear dynamics. Bulletin of Mathematical Biology, 85(11):110, 2023.
Johannes Brandstetter, Max Welling, and Daniel E Worrall. Lie point symmetry data augmentation for neural pde solvers. In International Conference on Machine Learning, pp. 2241 2256. PMLR, 2022.
Steven L. Brunton, Joshua I. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113 (15):3932-300) 3937, 2016. doi: 10.1073/pnas 15173......1111.....
Kathleen Champion, Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45):22445 22451, 2019. doi: 10.1073, 20.1073/pnss... 190699906999.16.
Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Solving and learning nonlinear pdes with gaussian processes. Journal of Computational Physics, 447:110668, 2021.
13

---

Published in Transactions on Machine Learning Research (03/2026)
Miles Cranmer. Interpretable machine learning for science with pysr and symbollcregrressoon..t 2023.
Miles Cranmer, Alvaro Sanchez Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Newral Information Processing Systems, volume 33,
pp. 17429-17442 Curran Associates, Inc., 2020.
Miles D. Cranmer, Rui Xu, Peter Battaglia, and Shirley Ho. Learning symbolic physics with graph networks, 2019.
David Dalton, Dirk Husmeier, and Hao Gao. Physics and lie symmetry informed gaussian processes. In Forty-first International Conference on Machine Learning, 2024.
Brian de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, Kutz, and Steven Brunton. Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49):2104, 2020. doi: 10.21 00.21105/joss.02100 URL https://doi.orr org/10 10. 21105/
jo83. 02104
Felix Dietrich, Alexei Makeev, George Kevrekidis, Nikolaos Evangelou, Tom Bertalan, Sebastian Reich, and Ioannis G Kevrekidis. Learning effective stochastic differential equations from microscopic simulations: Linking stochastic numerics to deep learning. Chaos: An Intentisciplinary Journal of Nonlinear Science, 33(2), 2023.
Renáta Duběáková. Eureqa: software review, 2011.
Gianluca Fabiani, Nikolaos Evangelou, Tianqi Cui, Juan M Bello-Rivas, Cristina P Martin-Linares, Con- stantinos Siettos, and Ioannis G Kevrekidis. Task-oriented machine learning surrogates for tipping points of baset models. Nature communications, 15(1):4117, 2024.
Mare Finzi, Gregory Benton, and Andrew G Wilson. Residual pathway priors for soft equivariance constraints. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 30037 30049. Curran Associates, Inc., 2021a. URL httpps://proceeddigg www.snnuuriippoooooeo pdf
Mare Finzi, Max Welling, and Andrew Gordon Wilson. A practical method for constructing equivariant multilayer perceptrons for arbitrary matrix groups. In International Conference Machine Learning, pp. 3318 3328. PMLR, 2021b.
Daniel Floryan and Michael D Graham. Data-driven discovery of intrinsic dynamics. Nalure Machine Intelligence, 4(12):1113 1120, 2022.
Evangelos Galaris, Gianluca Fabiani, Ioannis Gallos, Ioannis Kevrekidis, and Constantinos Siettos. Numerical bifurcation analysis of pdes from lattice boltzmann model simulations: parsimonious machine learning approach. Journal of Scientific Compuling, 92(2):34, 2022.
Sébastien Gaucel, Maarten Keijzer, Evelyne Lutton, and Alberto Tonda. Learning dynamical systems using standard symbolic regression. In Miguel Nicolau, Krzysztof Krawiec, Malcolm I. Heywood, Mauro Castelli, Pablo Carcía-Sánchez, Juan Merclo, Victor M. Rivas Santos, and Kevin Sim (ods.), Genetic Programming, pp. 25 36, Berlin, Heidelberg, 2014. Springer Berlin Heidelberg.
Jason Graham, K Kanov, XIA Yang, M Lee, N Malaya, CC Lalescu, R Burns, G Eyink, A Szalay, RD Moser, et al. A web services accessible database of turbulent channel flow and its use for testing a new integral wall model for les. Journal of Turbulence, 17(2):181 215, 2016.
Arya Grayeli, Atharva Sehgal, Omar Costilla Reyes, Miles Cranmer, and Swarat Chaudhuri. Symbolic regression with a learned concept library. Advances in Neural Information Processing Systems, 37:44678- 44709, 2024.
14

---

Published in Transactions on Machine Learning Research (03/2026)
Arthur Grundner, Tom Beucler, Pierre Gentine, and Veronika Eyring. Data-driven equation discovery of cloud cover parameterization. arXiv preprint or arXii:2304.00000 08063, 2023.
Daniel R Gurevich, Matthew R Golden, Patrick AK Reinbold, and Roman O Grigoriev. Learning fluid physics from highly turbulent data using sparse physics-informeet discovery of empirical relations (spider). Journal of Fluid Mechamics, 996:A25, 2024.
Samuel Holt, Zhaozhi Qian, and Mihaela der Schaar. Deep generative symbolic regression. Xre preprint arXiv:2401.00282, 2023.
Krzysztof Kacprzyk, Zhaozhi Qian, and Mihaela van der Schaar. D-cipher: discovery of closed-form partial differential equations. Advances in Newral Information Processing Systems, 36:27609-27644, 2023.
Kadierdan Kaheman, J Nathan Kutz, and Steven Brunton. Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A, 476(2242))20200279, 2020.
Pierre-Alexandre Kamicnny, Stéphane d' d'Ascoli, Guillaume Lample, and François Charton. End-to-end symbolic regression with translormers. Advances in Newral Information Processing Systems, 35:10269 10281, 2022.
Alan A. Kaptanoglu, Brian M. de Silva, Urban Fasel, Kadierdan Kaheman, Andy J. Goldschmidt, Jared Callaham, Charles B. Delahunt, Zachary G. Nicolaou, Kathleen Champion, Jean-Christophe Loiseau, J. Nathan Kutz, and Steven L. Brunton. Pysindy: A comprehensive python package for robust sparse system identification. Journal of Open Source Software, 7(69))3994, 2022. doi: 10.21105/joss.03996 03994. URL
https://doi.. www.ooor///////
12///10.22 2tt///00211155 21..2115//oss.. 03994
Gyeonghoon Ko, Hyunsu Kim, and Juho Lee. Learning infinitesimal generators of continuous symmetries from data. Advances in Neural Information Processing Systems, 37:85973 86003, 2024.
Kookjin Lee, Nathaniel Trask, and Panos Stinis. Structure-preserving sparse identification of nonlinear dynamics for data-driven modeling. In Bin Dong, Qianxiao Li, Lei Wang, and Zhi-Qin John Xu (eds.), Proceedings of Mathematical and Scientific Machine Learning, volume 190 of Proceedings of Machine Learning Research, pp. 65-80, PMLR, 15-17 Aug 2022.
Seungjoon Lee, Mahdi Kooshkbaghi, Konstantinos Spiliotis, Constantinos 1 Siettos, and loannis G Kevrekidis. Coarse-scale pdes from fine-scale observations via machine learning. Chaos: An Interdisciplinary Journal of Nontinear Science, 30(1), 2020.
Seungjoon Lee, Yorgos M Psarellis, Constantinos Siettos, and Ioannis G Kevrekidis. Learning black-and gray-box chemotactic pdes/closures from agent based monte carlo simulation data. Journal of Mathemationl Biology, 87(1):15, 2023.
Yi Li, Eric Perlman, Minping Wan, Yunke Yang, Charles Meneveau, Randal Burns, Shiyi Chen, Alexander Szalay, and Gregory Eyink. A public turbulence database cluster and applications to study lagrangian evolution of velocity increments in turbulence. Journal of Turbulence, (9):N31, 2008.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigode Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3():2118 229, 2021.
Georg Martius and Christoph Н Lampert. Extrapolation and learning equations. arXiv preprint arXiv:1610.02995 2016.
Mattoo Merler, Katsiaryna Haitsiukevich, Nicola Daineso, and Pekka Marttinen. In-context symbolic regression: Leveraging large language models for function discovery. arXiv preprint arXiv:2404.19094 2024.
15

---

Published in Transactions on Machine Learning Research (03/2026)
Daniel A Messenger and David M Bortz. Weak sindy for partial differential equations. Journal of Computa- lional Physics, 443:1 10525,, 2021a.
Daniel A Messenger and David M Bortz. Weak sindy: Galerkin-based data-driven model selection. Multiscale Modeling & Simulation, 19(3)14474 1497, 2021b.
Daniel A Messenger, Joshua W Burby, and David M Bortz. Coarse-graining hamiltonian systems using wsindy. Scientific Reports, 14(1):14457, 2024.
Grégoire Mialon, Quentin Garrido, Hannah Lawrence, Danyal Rehman, Yann LeCun, and Bobak Kiani. Self- supervised learning with lie symmetries for partial differential equations. Advances in Neural Information Processing Systems, 36:28973 29004, 2023.
Alan C Newell. Solitons in mathematics and plysics. SIAM, 1985.
Peter J Olver. Applications of Lie groups to differential equations, volume 107. Springer Science & Business Media, 1993.
Peter J Olver. Equivalence, invariants and symmetry. Cambridge University Press, 1995.
Samuel E. Otto, Nicholas Zolman, J. Nathan Kutz, and Steven L. Brunton. A unified framework to enforce, discover, and promote symmetry in machine learning, 2023.
Brenden K Petersen, Mikel Landajuela, T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne T Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXiv preprint orXiv:1912.0487. 2019.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihacla van der Schaar. D-code: observed trajectories. In International Conference on Learning Representations, 2022.
Discovering closed-form odes from
Chengping Rao, Pu Ren, Yang Liu, and Hao Sun. Discovering nonlinear pdes from scarce data with physics-cncoded learning. aixiv preprint arXiv:2201 12354. 2022.
Samuel Н Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz, Data-driven discovery partial differential equations. Science advances, ((4)=e6602614 2017.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. In International Conference Machine Learning, pp. 4442 4450. PMLR, 2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324 (5923):81-88) 2009.
Parshin Shojace, Kazem Mcidani, Shashank Gupta, Amir Barati Farimani, and Chandan K Reddy. Llm-sr: Scientific equation discovery via programming with large language models. arXie preprint (044:2400.18400, 2024.
Parshin Shojaee, Ngoe-Hieu Nguyen, Kazem Meidani, Amir Barati Farimani, Khoa D Doan, and Chandan K Reddy. Lim-srbench: new benchmark for scientific equation discovery with large langnage models. arXiv preprint Xiv::22044. 10415, 2025.
Makolo Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35:1596 1611, 2022.
April Tran and David Bortz. Weak form scientific machine learning: Test function construction for system identification. arXiv preprint 2XXi:2557.03200 2025.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16):eaay26)) 2020.
16

---

Published in Transactions on Machine Learning Research (03/2026)
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Newral Information Processing Systems, 33:4860-4811 2020.
Pantelis R Vlachas, Georgios Arampatzis, Caroline Uhler, and Petros Koumoutsakos. Multiscale simulations of complex systems by learning their effective dynamics. Noture Machine Intelligence, 4(4):359- 3666, 2022.
Pantelis-Rafail Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos, Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spaliolemporal dynamics. Neural Networks, 126:191- 217, 2020.
Rui Wang, Robin Walters, and Rose Yu. Incorporating symmelry into deep dynamics models for improved generalization. In International Conference on Learning Representations, 2021.
Rui Wang, Robin Walters, and Rose Yu. Approximately equivariant networks for imperfectly symmetric dynamics. In International Conference on Machine Learning. PMLR, 2022.
Yiqun Wang, Nicholas Wagner, and James M Rondinelli. Symbolic regression in materials science. MRS Communications, 9(3):793-8D0 2019.
Xiaoyu Xie, Arash Samaei, Jiachen Guo, Wing Kam Liu, and Zhengtao Gan. Data-driven discovery of dimensionless numbers and governing laws from scarce measurements. Nature Commenicatiinss, 13(1):5662 2022. doi: 10.1038/s4144 01/////4447-022-3555 22223508--..
Jianke Yang, Wang Rao, Nima Dehiaay,, Robin Walters, and Rose Yu. Symmetry-informed governing equation discovery. In Adconces in Neural Information Processing Systems (NeurIPS), 2024.
Zhi- Yong Zhang, Hui Zhang, Li-Sheng Zhang, and Lei-Lei Guo. Enforcing continuous symmetries in physics- informed neural network for solving forward and inverse problems of partial differential equations. Journal of Computational Physics, 492:112415, 2023.
17

---

Published in Transactions on Machine Learning Research (03/2026)
A Math
A.1 Notations
Table 2: Descriptions of symbols used throughout the paper. The three blocks include (1) basic notations for PDEs, (2) notations for Lie symmetry of PDEs, and (3) notations for symbolic regression algorithms and miscellancous.
Symbols
Descriptions
p
Number of independent variables of PDE.
q
Number of dependent variables of PDE.
X
Space of independent variables of PDE: X C IP. Also used to denote the feature space of SR algorithms.
U
Space of dependent variables of PDE: U C R Assumed to be 1-dimensional unless otherwise stated.
E
Total space of all variables of a PDE: E x X
Uk
Space of strictly kth-order partial derivatives of variables in U w.r.t variables X.
U
Space of all partial derivatives up to nth order (including the original variables in U): U(m) = U1 X Un-
((
ath-order jet space: M(a) c U(m).
TM
The tangent bundle of manifold M. Independent variables of PDE: E R
x
Time variable.
r,y
Spatial variables in PDE contexts. Also used to denote the features and labels of SR algorithms, where æ can denote multi-dimensional features.
u, u
Dependent variable(s) of PDE: u E and u E R?.
(m)) m
The collection all up to n-th order partial derivatives of or u. The (ordinary) differential of a function. For a differential function f ((n)) + RR df
df
En of 827 dx + E., Bua 0f dud
D.f
The total derivative of differential function f M(n) + w.r.t the ith independent af variable. For example, if p g = 1, Dif = äz + Σx= Us+180g where Ug == 8h u/ u/dx".
as
Dj
The total differential of a differential function f M() + R, i.c. Df = D,f dx'
g
A group element with an action on E (2).
v
vector field on the total space 3). representing an infinitesimal transformation. list of multiple vector fields are indexed by subscripts.
pr(n) pri"iv glv), v(n)
ath-order prolongation of g acting on M())
nth-order prolongation of v acting on M(a)
Equivalent to pr("), y and pr("),, v, respectively.
The (infinite) prolongation of v. For an nth-order differential function f(x, u(x)), pr v(f)- pf())
pr v
7): C, 8
Differential invariants of symmetry group. Y) is used by default. The other letters are used distinguish between invariants of different orders.
t,
The LHS of SINDy equation ( Often assumed to be time derivatives. A column vector containing all SINDy library functions: = [0] gm3
0
w, W
The SINDy parameters. For only one equation, w = [w2 ... w" is a row vector. For multiple equations, W == [w
9 x m matrix.
X,y
Concatenated matrix/ matrix/vectar /vctor of features /labels of all datapoints for symbolic regression.
[N]
List of positive integers up to N, i.c. [1,2,... ,]] for any N z
I:N
Equivalent to [N].
LHS, RHS
Left- and Right-hand side of an equation.
18

---

Published in Transactions on Machine Learning Research (03/2026)
A.2 Extended Background on PDE Symmetry
References for the below material include Olver
1993
Olver 1995
Prolonged group actions Let = x UUX = R x R be endowed with the action of group G via point transformations Then group elements E G act locally on functions u = [(x), therefore also on derivatives these functions. This in turn induces, at least pointwise, "prolonged" transformations on jet spaces: (x, (()) = pr(")g (x, (())..
Let J= (j,, -. du S p be n-tuple of indices of independent variables and 1 S a S q. We will use the shorthand
0'
ulu Oz,, 'z0
dx
and
DJ = Dj Dja-
It not practical to work explicitly with prolonged group transformations. Therefore one linearizes and considers the prolonged action of the infinitesimal generators of G. Explicitly, given a vector field
a Dua
Σ
(x,
+ Σ
(x, u)
Oz
a=1
its characteristic is a q-tuple Q (Q', ...., Q?) of functions with
(x (x, > Ex u) Oua
Ort
Now the prolongation of v to order is defined by
(n) = Σ e u)
Σ
Σ
(x,
u (i)))
(9)
dr
Oug
#Jmn
Here J ranges over all n-tuples == U......;n))
ju S p and the fo given by
+9 = D2Qa= + Σ e J.
We remark that the prolongation of v has been described explicitly in terms of the coefficients of v and their derivatives.
Figure 5 visualizes the group action of a Lie point transformation and its prolongation with a simple example. Consider the total space U = R, and the standard rotation generator in 2D space given by v= -ud + x0,, The vector field visualized in dark red arrows in the background. We also consider a function X U given by u(x) 0.5(x 1)3, whose graph is visualized by the dark red solid line in Figure 5 left. The graph of its first-order derivative, (x)) = 1.5(x 1)2, is visualized by the dark green dash-dot line.
Then, we choose a random group element g = exp(Ov) that rotates a 2D vector (x, u) E X x U by angle 0. Applying this pointwise transformation to every point on the graph of u(x), we have a transformed graph visualized by the dark red dashed line. The transformed function, = g u, is defined as the function whose graph is the transformed graph. In other words, ū(x) = (g (x)x) visualized by the dark red dashod line Figure left.
du The prolonged Next, consider how the rotation of (z,u) transforms the first-order derivative u de vector field, ie., the infinitesimal generator of the prolonged group action, can be computed by ( v
+ (1 + uZ) 0 The projection of v(1) onto X Un is visualized in the dark green arrows in Figure
19

---

Published in Transactions on Machine Learning Research (03/2026)
Rotation (g) of u(x) and its prolongation to Uy(x)
Full view of g in X x U x Ux =
Durves
Original utzl,
uix)
Transformad ed curve: (E, 900) 0,09)
Gix) (g- -100
Original pro ected
uyD0
Transformed ed projected to X x and X
(xlx) -u.ld)
Vi obyy XOu
1/2 -B, +
Ux
x
Figure 5: Demonstration of the rotation v -ud, xún acting on x U, and first-order prolongation acting on U U1.
left. Similarly, the prolonged group action gL2 = (x)(Ov(1)) is applied to every point on the graph of un((), dē yielding the graph of the transformed derivative function, üz (x) du (z), visualized in the dark green dotted line.
The full transformation ion of the prolonged (4) in the 3D space X U x Un is shown Figure right. The graph of the original prolonged function (1)) (a)) (u(x)), (x())) is shown in the solid line, which transformed into the dashed line by (2))
A.3 Proof of Proposition
4.3
Olver 1995 provides the following general theorem to construct higher-order differential invariants from contact-invarianl coframe. We refer the readers to Chapter 5 Olver (1995) for definitions of relevant. concepts, e.g., contact forms and contact-invariant forms and coframes.
Theorem A.1 (Thm. 5.48, (Olver 1995)). Let G be transformation group acting on a space with independent variables and q dependent variables. Suppose w',...,'' wD is a contact-invariant coframe for G, and Dj the associaled invariand differential operators defined via DJ = D2J dx = =D, w. If there are a number of independent, strictly nthh-order differential invoriants <².... 5"h, qn = (p-n-1). then the set of differentiated invoriants D,S", i E [p), [an], contains complete set of independent, strictly (n+ 1)th-order differential invariants.
Specifically, the condition that there exist maximal number of differential invariants of order exactly is guaranteed if n is at least dimG.
Our proposition derived result from the above theorem, which provides a concrete way of computation from lower-order invariants higher-order ones:
Proposition A.2. Let G be local group acting on X R2 R. Let 1/², "..."..,,,. be any p differential invoriants of G whose horizontal Jacobion J [Dop'] is nen-degeneaate on an open subset 2 C (") there are marimal number of independent, strictly nth-order differential invariants ç!.... _6 On (P+n-1). then the following sel contains complete set of independent, striclly (n + 1)th-order differential invariands defined on Ω:
n
dct(D. Wk E [p], k' E [gn),
(10)
det det(Den²)
20

---

Published in Transactions on Machine Learning Research (03/2026)
where i,j E [p] are matrir indices, Dj denotes the total derivative w.r.t i-th independent variable and ((B,))) = "...,,,-1,, , 1.... "/.+,.....?
Proof. We show that the total differentials of the differential invariants ry',... 7p can be used to construct a contact-invarint coframe G and then derive the associated invariant differential operators to complete the proof.
First, note that for any differential invariant y of G, its total differential Dn Dyn dx can be written
we + 0,
(11)
where w, dr- Cie(p) de dx Ejaggn
OF
OF dua is the ordinary differential of r M(x) IR and 0 is Du.
contact form.
Since is differential invariant, its differential w. dn is an invariant one-form on M M(),, i.e. (g wor
Also, a prolonged group action maps contact forms to contact forms. To see this, note that a prolonged group action gln) maps the prolonged graph of any function to the prolonged graph of a transformed function. Then, for any contact form 8, (g("))" annihilated by all prolonged functions f(n), thus a contact form by definition:
(( (g(") piuly00
= ((g (I (1))00
(12)
Then, from 11
we have
(g(atl))" = )" + (g (1
wo + f
(0' 0)
(13)
where 8' is some contact form and so is 0' 0. Thus, is contact-invariant. For the p differential invariants .., ,77,, have p inntactt-nvvaiant one-forms w/, .w, respectivcly.
Next, we prove that w'.... .... w are linearly independent and form coframe. Assume there exists smooth coefficients e² such that E, we = 0. Then, regrouping the coefficients of the horizontal forms dx', we have
0 Σ Dam 1
Σ
Σ
è Dar de'.
(14)
Because the dx are linearly independent, each coefficient. of do must vanish, i.e. J = Since the
Jacobian = [D,r'] is non-degenerate, the only solution is ~ = (on the open subset St E (m)).. Thus, w!.... " form a contact-invariant coframe. According to Thoorem A.1 the associated invariant differential operators of the coframe take complete set of same-order invariants invariants.
complete set of cae-order--ighee
The remaining step is to obtain the invariant differential operators explicitly in terms of n². Recall the formula in Theorem A. that defines the invariant differential operators:
Duf dx D,f wt
(15)
Expanding w = Dip = Dap dx', we have the following linear system of invariant differential operators Dy;
Di
Din'
Din²
DiT D21
D
D2
Dzn' Dzn²
D2
(16)
:
Dp
Dant
Dan²
Dam
Dp
21

---

Published in Transactions on Machine Learning Research (03/2026)
Since J [Dan'] is non-degenerate, Cramer's rule yields
D&ç= det(Dan)
Dant-1
| Di | ... Danp) det(DDp))
(17)
Remark A.3. We require that the differential invariants ....... ,7 has nondegenerate horizontal Jacobian [Drr',, which is stronger condition than functional independence. Since the differential invariants are functions on the jet space, it is possible that a set of such functions is functionally independent, i.e, has a nondegenerate full Jacobian [0,r'], where E [qn] indexes the jet space variables (x, ), but has a lower-rank horizontal Jacobian. For example, consider ry' u, and 1y² Ny- In the full Jacobian, Önp /0u, and 0np' /Oug form the identity, it has full rank. However, its horizontal Jacobian containing total derivatives given by Maz
Vay
which is not invertible on the subset of the jet space where Unr Uzy
" 0.
May
In practice, this non-degeneracy condition can be easily checkod once we have the symbolic expressions of the p differential invariants.
Remark A.4. When 1, Proposition A.2 is equivalent to the following (Prop. 2.53, Olver 1993)
If y = 7x,, ( and w = Ç(x, ( are n-th order differential invariants G, then du = D.S is an (n + 1)- th order differential invariant of G. Specifically, if = 77, u) and w <(x,, u, Uz) form complete set of functionally independent differential invariants of pr the complete set of functionally independent differential invariants for pr is then given by
Dz7
y, w. dw/ dy, ... d= w/ /dy"d""
(18)
A.4 Examples of Computing Differential Invariants
Example A.5. Consider the group SO(2) acting on X XxU=R² x R by standard rotation in the 2D space of independent variable and trivial action on U, Le. its infinitesimal generator given by v= yy2 - xd,.
First, we solve for complete invariants 7 1), y, u) should satisfy yd,7 z0,7 = 0. Since the vector field does not involve u, immediate solution is u. On the othe hand, by method of characteristics, we convert the PDE to the characteristic equations dr /r/ds= y, dy/ds -r. That is, the characteristics curves (I(s), ((s) are just circles around origin. Because is constant along characteristic curves, it must be a function of R²2 = x² + +1² Therefore, we pick the following two ordinary invariants: (1, y, u) = [x² + y²) and (2 y. u) u. (5) dictates how we construct higher-order invariants using these two functionally independent invariants and another arbitrary invariant. For notational convenience, we convert to operators defined according 7/2 and m respectively:
of the ordinary and first-order invariants. By definition, the ordinary
O xDyy yDz
(19)
Flly yu
O₂ Uy Dr Ug Dv
(20)
Fuy yu
Then, need to find another new differential invariant, because applying these operators on m and 7/2 leads to trivial results. Since 77 and 7/ generate all ordinary (zeroth-order) invariants, we must look for the first-order invariants. To do this, note the prolonged vector field is given by
pr ( v + uy - -uz
(21)
Solving for pr gives two first-order invariants, Çi vuy yuz and Ço xus + yug. Note that the differential invariant Çi is exactly the common denominator in O and O2, we can simplify O, and O₂ by
22

---

Published in Transactions on Machine Learning Research (03/2026)
using only their numerators, i.c.
0 = xDy - yD:
(22)
O₂ ug D= U2 Dy
(23)
Note that O2 has lirst-order coefficients, which may complicate things in the subsequent calculation. Denoting the space of all contimuus functions the existing four invariants as I C(n, m2, Ç1, (2), we can choose any new operator within the I-module spanned by O and O₂ that makes things casier. Specifically, we use the following operator
O 42 00 2th O2
(1 21
= xD + yDv
(24)
Then, we apply these operators the first-order invariants, which raise the order by one and give us the second-order invariants. For example, applying O Si, we have
O11I = IDySI - yD
= xxuuyy - ur yuzu) - y(u, + Iuly - yllz)
= unu - y² uex Tul yuly 2ryuay
(25)
Note that (2 vu, + yiyy is first-order invariant, so can further remove from the formula and get simplified second-order invariant
91 = x² uyo + y Unn 2xyuuy
(26)
Similarly, we compute O1(2, O2h1 and Ö2L2 and obtain the following, respectively:
022 03 S1 + zy Uxx) (x² - yy2u
= ryyluv Uxx) (x² 15² yuuu
(27)
04 = C2 + x² x²ugx + uu + 2xzyuz
= x uge + yuugg + 2zyuxy
(28)
The above invariants should form complete set of second-order differential invariants of v x. yd To verify, note that the Laplacian Au= Uez + Uyys which is well-known rotational invariant, can be written in terms of these differential functions:
(x² + y/2) (Uaa + Ugg
Au Ugz + Uyy
x² + yz²
01 04
(29)
21/1
Another second-order rotational invariant, the trace of the squared Hessian matrix, u2 + 2u + u2 recovered by
wy: is
0 07 20 + 02
22 + 2u2 + An
(30)
4m
On the other hand, these 8 invariants are apparently not functionally independent - note that 12 = O and 0g ŌzCi are the same. While this may some coincidence, eventually it is not surprising because we would expect see 3 functionally independent strictly second-order differential invariants instead of 4, since (Uzz, Myg. Mny) U2 is only 3-dimensional.
23

---

Published in Transactions on Machine Learning Research (03/2026)
Example A.6 (Scaling and translation). Consider the vector field V1 tô, + axd, + buQ,. It generates the scaling symmetry 1 + Mr >> x² + A'u. The ordinary invariants of this symmetry are t'u and w0=-1. The higher-order invariants are given by [a,8) = u where and 8 denote the orders of partial derivatives w.r.t t and z, c.g. "2(ayy0 == Uazt-
Besides the scaling symmetry, we can consider other common symmetries simultaneously, e-g. translation symmetries in both space and time, V2 =2 = 0g and V3 = These symmetries, along with the scaling symmetry VI, span three-dimensional symmetry group. There are no ordinary invariants due to the translation symmetries. A convenient maximal set of functionally independent differential invariants is given by
(a.3) = "g
a 2 0,
3 2 0.
(31)
A.5 Proof of Proposition 4.4
Proposition 4.4 restated below, aligus our symmetry constraint into the SINDy framework and results in a set of constraints on the SINDy parameters.
Proposition A.7. Let l(x, (*)) = W u(*)) be a system of differential equations admitting a symeeerr group G, where x E R2, u E R2, E IRm Assume there erist some nih-order invariants of G, Yl 1:9 and 17" n2 s.l. (1) the system of equations can be expressed as no = W')')n) where vp- = [1o 1:9] and n = [,j=*], and (2) no = T and (0')² = SS) for some functions (n) and constant tensors W and Then, the space of all possible W is linear subspace of Roxm.
Proof. (Note: In this proof, we do not distinguish between superscripts and subscripts. All are used for tensor indices, not partial derivatives.)
For simplicity, we omit the dependency of functions and write
= Wr
(32)
Combining the conditions about the differential invariants, know that the equation can be equivalently expressed as
ghi ei (W')'sgske)
(33)
for some W' E R2xm', where m' is the number of invariant functions in 0'.
32 into 33 and rearranging the indices, the principle of symmetry invariants then translates to the following constraint on W: there exists some W E Roxm' s.t.
Substituting
Tr 00 = (W') s Vx, (m).
(34)
To solve for W, we first eliminate the dependency on the variables and (n) from the equation. We adopt procedure similar to Yang et al. (2024). Denote z (x, (*)). Define functional Me mapping function to its coordinate in the function space spanned by 0, i.e. Me (z + cz0) + (c),c² ....c...... ..").. Before we proceed, note that the LHS of (34) contains the products of functions (z)) which may or may not be included in the original function library
Therefore, we denote 0(2) = (0(z) || (0x0: & 0)] as the collection of all library functions 0₁ and all their products 0g0,. The invariant functions 0' (((1)) can also be rewrilten terms of the prolonged library: 0' ('(n) S0, where S1. S.
Then, applying Mg
34
have
(.TT, = 3.
(35)
Further expanding the LHS, we have
T W =(W==')''3
(36)
24

---

Published in Transactions on Machine Learning Research (03/2026)
where I satisfies 0n0: = 00). In other words, the rows of the LHS fall in the space of S. Let S1 be the basis matrix for the null space of S, i... 85- = O, we have
T w, " ( ) 0, suggesting that W must lie in a linear subspace of R2xm.
(37)
we first rearrange (37 into Mv) = O, where M has shape (S.shape(2) x 9, q m). Then, we perform SVD on M and apply threshold of 10-6 6 to the singular values. The right singular vectors corresponding to the singular values smaller than the threshold then form basis of the linear subspace vec(W) lies in.
Remark A.8. In practice, to solve for 37
B Implementation Details
This section discusses some detailed considerations in implementing the sparse regression-based methods described in Section 4.3 and
Contents include:
Appendix B.1
An algorithmic description of direct sparse regression with symmetry invariants.
Appendix B.2
Converting the symmetry invariant condition as linear constraints on the sparse
regression parameters.
Appendix B.3 Using differential invariants in weak SINDy via the linear constraints, as well as other
considerations.
B.1 Direct Sparse Regression With Symmetry Invariants
The first approach enforcing symmetry in sparse regression, as discussed in Scction 4.3 is directly use the symmetry invariants the variables and their functions specified by a function library the RHS features. Similar to Algorithm | for general symbolic regression methods, we provide detailed algorithm for sparse regression below. Following the setup from SINDy, we aim to discover system of g differential equations for q dependent variables.
Algorithm Sparse regression with symmetry invariants
Require: PDE order datasct fz' = (x', (o(),((())))) (u(*))) E M( li=1: SINDy LHS e, SINDy function library (0),
infinitesimal generators of the symmetry group B = (va)-
Ensure: A PDE system admitting the given symmetry group.
Compute the symmctry invariants of B up nth-order: ,,.,... k (Prop. 4.3
Choose an invariant function 1² s.t. Ory²: /e * 0 for SINDY LHS component
Tet no [,)x,, ",.,,,..,,> "T and v denote the column vector containing the remaining K - invariants.
Instantiate the sparse regression model 770 W'()).
+| return 70 = We(n). [Optionally, expand all n in termss of original variables z.}
Optimize W with the SINDy objective: Σ: [1') ) -
The configuration from the original SINDy model, i.e., the LHS & and the function library {0) }, are used to construct a new equation model in terms of the invariants. should be noted that the functions in the SINDy function library does not specify their input variables. For example, in the PySINDY (Kaptanoglu et al., b022) implementation, a function provided a lambda format lambda x, y: x y. Thus, O can be applied to both the original variables, e.g. 0(21,2) z1 and the invariant functions, e.g. 0(n, 1//2.
B.2 Symmetry Invariant Condition Linear Constraints
Instead of directly using the invariant functions n as the features and labels for regression, we can derive a set of linear constraints from the fact that the equation be rewritten in terms of invariant functions. As
25

---

Published in Transactions on Machine Learning Research (03/2026)
shown in Appendix A.5 basis Q of the constrained parameter space can be obtained from the right singular vectors of a constraint matrix M. We rearrange Q to a tensor of shape (r, m)), where r is the dimension of the constrained parameter space, and (q, m) is the original shape of the parameter matrix W. Then, we can parameterize W by Wik = Q' 3t, where / is the learnable parameter, and discover the equation using the original SINDy objective as described in Section 3.2
In practice, we observe that the basis Q obtained from SVD on 37 is not sparse. Indeed, SVD does not inherently encourage sparsity in the singular vectors. As mentioned in Section [.3] the lack of sparsity pose a problem when we perform sequential thresholding in sparse regression. Therefore, after performing SVD, we apply a Sparse PCA to Q to obtain a sparsified basis, also of shape (r; (r.g. m))
spca =
(n_ commooonttt_)
spca. fit (Q. reshape q*m)) Q_ sparse spca. components reshape
Figure Q shows an example of the original basis solved from SVD (top 7 x grid) and the sparsified basis using sparse PCA (bottom 7 x grid). This is used in our experiment the reaction-difusion system (8)
B.3 Using Differential Invariants in Weak SINDY
In this subscction, discuss the formulation of weak SINDy and how to implement our strategy of using differential invariants within the weak SINDy framework. To maintain similar notation to the original works on weak SINDy (Messenger & Bortz 2021a we use Do, to denote partial derivative operators, where a₂ = (81, $2...... Sp) is a multi-index, instead of using subscripts for partial derivatives. Thus, we no longer strictly differentiate subscripts and superscripts- botth can be used for indexing lists, vectors, etc.
Given a differential equation in the form
Dao u
Σ Wa.sDa.lJuu 8.
(38)
can perform integration by parts (i.c., divergence theorem) to move the derivatives from to some analytic test function and thus bypass the need to estimate derivatives numerically. First, multiply both sides of 38 by test function o with compact support B C. X and integrate over the spacetime domain:
J Dan u(x)x(x))xx
Σwa, Waj / D (xxx))
(39)
s.j
WLOG, assume that 81 * O, and denote as' = (s1 - s2,, Then, each term in the RHS can be integrated by parts as
J Da. (ux))((())))(
/ Da fy (u(x))xxx))
Da, (al,)((ul))))))(xx)d / in
=
[ Dul,f(u
(40)
where Dj denotes the partial derivative operator w.r.t the first. independent variable, and in is the first component of the unit outward normal vector.
Repeating this process until all the derivative operations move from fy( the test function o, we have
J Da (xxx (--)))/x)xx(())o))))) / ((,
(41)
Similarly for the LHS:
[ Dan ((xx)xxxxxxx(xII)))) / ((x))an((x))))
(42)
26

---

Published in Transactions on Machine Learning Research (03/2026)
d: du tp
oit shasullenniin dt 67 Ju A.r at
abd (ntettonooonnneee +52 40.0 4.2 43 0.1 0L9 233 4333
oi diw.siieeneees ai de aur at
Figure 6: Basis for the SINDy parameter subspace that preserves S0(2) symmetry SINDy parameter W has dimension 2 19. The two correspond to the two equations with u, and as the LHSs. The RHS contains 19 features, including all monomials of up to degree 3 and their spatial derivatives up to order 2. The set of symmetry invariants used to compute the basis is given by [t, x. y, v² + x²} Ufu u Ufu un where == (u, u)" and a multinnex t, x. with order no more than The top 2 grid displays the original basis solved from SVD, and the bottom the sparsified basis.
-vd, + ud,, The
to
grid displays
27

---

Published in Transactions on Machine Learning Research (03/2026)
The final optimization problem is to solve for b Gw, where w is the vectorized coefficient matrix W, and each row in b and G is given by computing the integrals in 441p and (42) against. a single test function. The number of rows equals the number of different test functions used.
Direct integration of symmetry via linear constraints As we have discussed in Appendix B.2 we can enforce symmetry by converting it to set of linear constraints the parameter W. With this approach, we can directly incorporate symmetry in weak SINDy. Specifically, we just parameterize W as in terms of a precomputed basis Q and a trainable vector 3 and directly substitute this parameterization of W into the optimization problem of weak SINDy. We adopt this strategy in our experiments concerning weak SINDy.
Expressing the equations with differential invariants The above approach is only possible when the conditions in Proposition 4.4 about the selected set of symmetry invariants hold. We should note that it is not always possible to find set of invariants so that the symmetry condition can be converted to linear constraints on the parameter W via the procedure in the proof of Proposition 4.4 One may ask the following question: can we simply express the equations in termss of differential invariants and apply weak SINDy, similar to Algorithm 2 for the original SINDy formulation? Here, we do not provide definite conclusion for this question, but only discuss several cases where directly using differential invariants in equations might succeed or fail in weak SINDy.
To adapt to the weak SINDy formulation 38 it is more helpful to consider the symmmelry invariants generated by some fundamental invariants and some invariant differential operators, instead of specifying a complete set of differential invariants for every order. Concretely, there exists a set of invariant differential operators } and a set of fundamental differential invariants I = {mHx s.t. every differential invariant can be written as jj...Oj. The For the SO( symmetry group in Example A.5 one possible choice is
7 2 ==2 +2²), 7/2 u, O, = =x,, yDa, O2=xDz x0=xx,+DD,
(43)
We can compose these generaling invariant operators to obtain a full library of eligible differential operators up to some order, denoted D = =D The exact compositions can vary and we can choose the most convenient one for subsequent calculations. For the above SO(2) example, for up to socond-order differential operators, we can choose O2, 0f. 02, v (O² + 02)}. Note the last operator exactly the Laplacian.
Then, the complete set of eligible terms (respecting the symmetry) in the equation is {D,ms : DD E D. E IJ. If assume, as in SINDy, that the governing equation can be written in linear combination these symmetry invariants, then can assign weight for each D, and form coefficient matrix W [W,x]- That is,
Djo lko
Σ
WSK nk
(44)
(io,))
Then, multiplying each side by a test function ((x), we have
f Dio (xx)dxx
Wgk / Dynos((xddx.
Σ
(45)
(..)7)ookk))
The question then boils down to whether we can apply the technique of integration by parts similarly to this set of differential operators and differential functions, since the original algorithm only deals with partial derivative operators Da, and ordinary functions (u).
To check this, let explicitly write out the dependency of these operators and fundamental invariants.
28

---

Published in Transactions on Machine Learning Research (03/2026)
Case 1 A relatively simple case is when all invariant operators take the form D, >5=E (x) and 1/K ((xx). Each term in the RHS of (45) can be expanded
J Dyne(xx)dx
Σ / as (x(x) Mr Mz ((x))
- / (o.
(46)
Evaluating [16] does not require estimating partial derivatives of u. Therefore, weak SINDy can be applied to this case quite straightforwardly.
Case However, it is not always possible to have all Dj classical linear differential operators and all Th ordinary functions. For instance, in Example A.6 there are ordinary symmetry invariants due to the constraint of translation symmetry.
If we still have linear operators Dj = E. a. (x) but on the other hand we have differential functions Ti = mh (x, (l))), we can still perform integration by parts as (46), but the final result becomes
Σ(-)) / ne ul() (as(aa (xxxxxxxx)
(47)
meaning still have to evaluate whatever partial derivatives remain in M7-- It is possible that decrease the order of partial derivatives compared to vanilla sparse regression, but we cannot eliminate all partial derivatives compared to Weak SINDy without any symmetry information.
Case The most challenging case is when the invariant differential operators explicitly involve the partial derivative, such as Dj= = E, az(x, u(n) Doy. Then, similar to (47 integration by parts yields:
/ ne ne(x, ul [ax(x, (xxxxxx
(-1)lool
(48)
In this case, still need to compute the partial derivatives, not only those in Vik but also those arising from a, and Da, (ag). The latter might involve higher-order derivatives and the benefit of using the weak formulation may further diminish.
C Additional Experiment Results
Contents of this section include:
Appendix C.1 Extended results in Table 1 with confidence intervals for the prediction error metric
over 100 runs.
Appendix
Results for some variants of the sparse regression models considered in Table I
Appendix
Results for the symmetry-breaalee reaction-diffusion systems with auto-tuned WSINDy parameters, complementing Figure?
Appendix C.4
Results on a 3D PDE system with SO(3) spatial symmetry.
Appendix
Results for the D-CIPHER (Kacprzyk al. 2023 baseline and our method applied to D-CIPHER the Darcy flow dataset.
Appendix
C.6
Samples of equations discovered by different methods.
Appendix
C.7
Visualized prediction errors of equations discovered by different methods.
29

---

Published in Transactions on Machine Learning Research (03/2026)
C.1 Results in Table 1 With Error Estimates
Table 3:
Extended results in Table with confidence intervals for the prediction error metric over 100 runs.
Each table entry is formatted as median [25% quantile, 75% quantile].
Method
Boussinesq
Darcy flow
Reaction-diffusion
Sparse
PySINDy
0.373 [0.367, 0.380]
0.021 [0.020, 0.022
Regression
SI
0.098
[0.098, 0.098
0.008 [0.007, 0.013
Genetic
PySR
0.098
[0.098, 0.098
0.114 [0.089, 0.169
Programming
SI
0.098
(0.098, 0.098
0.051 [0.031, 0.053
0.023 [0.015, 0.036)
E2E
0.132
[0.109 0.322
Transformer
ST
0.104 [0.100, 0.109]
C.2
Variant Sparse Regression Models
Table 4: Results of sparse regression models the Boussinesq equation and the reaction-diffusion system. C stands for complexity, i.e, the dimensionality of the parameter space. SP stands for success probability. The PySINDy and SI rows present the same results as the corresponding rows in Table
Boussinesq
6
Reaction-diffusion
Method
C +
SP
C
SP
PySINDY PySINDy*
15
0.00
38
0.53
21
1.00
468
0.00
PySINDy **
15
1.00
198
0.00
SI
15
1.00
28
0.54
ST-aligned
14
0.56
PySINDy de Silva et al. 2020 Kaptanoglu et al. 2022) constructs its library O from list of variables and derivatives, [u] ua (la) > 0) and set of scalar functions specified in lambda format. For example, to specify the following functions: - and
include up to quadratic monomial terms in the library, (x, y) + zy. However, their original implementation does not allow these functions to be applied to partial derivative terms. As a result, terms such as u cannot be modcled. This leads to its failure to discover the Boussinesq equation (6) as we have shown Table
We modify the implementation and include an additional set of results with different libraries, denoted PySINDy* in Table 4 The PySINDy* model supports wider range of library functions, including functions of partial derivatives, e.g., uz. Further more, we notice that the PySINDy* library while comprehensive, contains many redundant terms, such as interactions between derivatives like ux Therefore, we implement another library, denoted PySINDy**, where functions such as (z, v) + xy are only applied when their arguments do not contain at least two different partial derivatives. Therefore, PySINDy"* library still includes all the necessary terms to recover the Boussinesq equation but becomes much more compact. A complete description of the hypothesis spaces of different sparse regression-based methods is available in Appendix E.5
As Table 4 shows, both PySINDy* and PySINDy** succeed in the Boussinesq equation. However, they fails in the reaction-diffusion system because their parameter spaces become too large due to a higher-dimensiooall total space X x U ~ R² x 2 Even with the more compact PySINDy**, there are still 198 possible terms for the reaction-diflusion system, and the algorithm never succeeded in 100 runs. This augments the point that SLNDy's success relies on an appropriate choice of function library. If the library is too small to contain all the terms appearing in the equation of interest, the discovery sure to fail. If the library too large, the optimization problem becomes more difficult the dinensimnniional parameter space. On the other hand, by introducing the inductive bias of symmetry, our method automatically identifies proper function library
30

---

Published in Transactions on Machine Learning Research (03/2026)
that contains all the necessary terms for PDE with specific symmetry group, but not other redundant terms.
We include another model in Table | Sl-aligned, where we derive a set of linear constraints on the sparse regression parameters from the fact that the equations can be expressed in terms of symmetry invariants. In this way, we still optimize the original parameters (though a constrained subspace) as the base SINDy model without symmetry, effectively "aligning" the hypotheses about equations from symmetry and the base SINDy model. This method is discussed in detail in Section 4.3 and Appendix [B.2] We should also note that this method is mainly developed for incorporating the symmetry constraints into the weak formulation of SINDy. However, it is perfectly acceptable to implement it in the original formulation of SINDy, so we provide its results in Table for reference.
For the reaction-diffssin system, Sl-aligned has a 14-dimensional parameter space. The basis for its parameter space is visualized in Figure 6 It achioves a slightly higher success probability than SI (regression with symmetry invariants) and PySINDy (without symmelry information). We do not apply SI-aligned to the Boussinesq equation, because it is not necessary to align the hypotheses from SINDy and symmetry in this case. We can readily convert any equation discovered from S1 (regression with symmetry invariants) by multiplying both sides by u
We note that the results on the reaction-diffusion system in Table are for models with the original SINDy formulation, in contrast to the weak SINDy formulation used in Figure Therefore, the results in Figure
should not be directly compared to those in Table and Table
C.3 Parameter Selection for Weak SINDY (WSINDy)
In Section 5.4 we have mentioned that different methods for selecting the test functions used in WSINDy have a significant impact on the equation discovery outcome. As is shown in Figure 4 (left), the auto-tuned test function parameters from the data result in better performance for baseline WSINDy (no symmetry constraint) in all noise levels considered, and for rrmmetry-constraine method under relatively low noise levels.
In Figure 4 (center & right), we have used the default test function parameters from PySINDy for the variants of the reaction-diffusion system with imperfect symmetry. We find that the default parameters result in higher success probabilities of finding the correct equations in these cases. However, for comparison, we also show the results of the same experiments using the auto-selected test function parameters in Figure
R-D w/ unequal diffusivities, autotuned WSINDY
R-D w/ external forcing, autotuned WSINDy
Si-relaxed
1.0
WSINDY
Perfect symmetry
Success probabiity
Success probability
0.6
0.6
Sl-relaxed
0.2
02
SI
WSINDY
Perfect symmetry
0.0
0.0
-0.03
-0.02
-0.01
0.01
0.03
Symmetry breaking parameter E
Symmetry breaking parameter E
Figure 7: Success probabilities on the cafftinndiifusion system with unoqual diffusivitics (left) and external forcing (right) with auto-tuned WSINDy parameters.
With the auto tuning of test function parameters following the procedure in Messenger & Bortz (2021a), the relaxed version of our method (SI-relaxed) still has superior performance the baseline WSLNDy in both
31

---

Published in Transactions on Machine Learning Research (03/2026)
symmetry- scenarios. Notably, the auto-tuning procedure typically suggests more test functions with higher polynomial degrees than the default setup, the performance of SI-relaxed is even better than its However, the baseline WSINDy does not benefit from this selection strategy
performance shown in Figurel ! as much, while the version of our method that enforces hard constraints fails to recover any equation with imperfect symmelry.
C.4 Equations with More Spatial Dimensions
Our previous examples focus on equations with 2 spatial dimensions at most. In this subsection, we test whether our method works well for discovering equations with higher spatial dimensionality. We consider a 3D isotropic reaction-diffusion system governed by
Ur = 0.2Au + 1I u.
(49)
where u is a scalar field, and the nonlinear reaction term is given by u u.
For data generation, we use a random perturbation around u = 0. Then, similar to the 2D reaction-diffusion system [8] we simulate the solution from this initial condition using an FFT-based spectral method with RK45 solver and periodic boundary conditions. We use 32 grid points per dimension and simulate up to T = 10, recording the solution at every A 0.1. After simulation, we add a 1% noise to the trajectory and compute the derivatives from the noisy trajectory. We then split the data along the time dimension, using E [0,8) for training and & [8, 10) for testing.
Table 5: Results on the 3D reaftion-diffusion dataset over 100 runs for each method. The prediction error (PE) column is formatted as median [25% quantile, 75% quantile].
Method
C 4
SP t
PE 1
Sparse
P PySINDy
286
0.51
0.216
[0.207, 0.977]
regression
S1
84
1.00
0.208
[0.204, 0.210
Genetic
PySR
10
0.00
1.228
[1.227, 1.230
programming
SI
6
0.20
1.227
[0.324, 1.229]
The equation admits SO(3) spatial symmetry. We use set of SO(3) invariants up to second order given by IRD-3D (u, (Vu)², Au, (Vu)² HVu, (r(HH),dot(H)) where H denotes the Hessian matrix HLy Uajs and ir and det stands for matrix trace and determinant, respectively. Note that a complete set of invariants would also include terms with independent variables and time derivatives, e.g., {ue, Uer ||x|| x-V,, x" Hx, x" Vu). However, the standard setup of SINDy excludes time derivatives and independent variables in the RHS features. To align with the SINDy setup, we only use those 6 invariants in IRD-ID. Similarly, for genetic programming, even though is possible to also include the spatial variables, we only use the standard coordinate in the 10-dimensional jet space (excluding t): {u, Uzi My: U2 Ugas Ugyr. Mya: Uzz). Finally, for the SINDy function library, we include all monomials up to degree 3 of all variables and derivatives.
Table 5 shows the results for sparse regression and GP-based methods on this 3D dataset. Our symmetry- seariant-based method (SI), when applied to SINDy, always discovers the correct functional form of the RHS of (49). In comparison, SINDy with regular variables only succeeds in about half of the 100 runs. Also, the complexity (C), which stands for the total number of feature terms for SINDy-based methods, is much lower with our method. As such, we comment that our method's advantage of reducing the equation space complexity becomes more obvious in these higher-dimensional examples.
On the other hand, genetic programming proves to be less favorable base SR algorithm in this case due to its overly large search space. PySR with standard jet variables completely fails to recover the correct equation. Our hod uses the invariant set TRD-3D as PySR input features and achieves a success probability (SP) of 20%. Finally, the symbolic transformer, which we used in the main experiments in Table I (E2E), has 0 success probability in this dataset, with either regular variables or symmetry invariants. We therefore do not report its results in Table 3 Again, we comment that the failure of the symbolic transformer might be because the data distribution from the PDE solution is largely different from its pretraining dataset.
32

---

Published in Transactions on Machine Learning Research (03/2026)
C.5 Comparison with D-CIPHER
A main advantage of our proposed method is its compatibility with various algorithms for symbolic regression. In the main experiments in Section 5.3 in Section 5.4 have shown that our method works well with SINDy (Brunton et al. 2016), weak SINDy (Messenger & Bortz 2021a), genetic programming (Cranmer 2023), and symbolic transformer (Kamienny et al., 2022). To further demonstrate this advantage, we include another base algorithm for symbolic regression, D-CIPHER (Kacprzyk et al. 2023), in this section.
Similar to weak SINDy, D-CIPHER (Kacprzyk et al.
2023 uses a variational objective for equation discovery.
It defines the extended derivative as
Elu](x)
u),
where a and h are some functions and a is a multi-index indicating partial derivatives. Then, a library {E") of such extended derivatives specified by the user by providing S triples of (a, a, h). The algorithm then optimizes for a coefficient / E R$ and a symbolic function g(x, u) under a variational loss and outputs the equation
Σ [u(x g(x, u).
As discussed in Appendix B33 our approach can be directly used in D-CIPHER to enforce symmetry. We demonstrate this with an experiment on the Darey flow dataset with S0( symmetry. First, we obtain the generating invariant operators of SO(2), i.c. 01 zD, yD, and O₂ xD, + yDy. To discover this second-order PDF, enumerate the following 5 up to second-order differential operators composed by O and O2: O = {O1, O2, Of, 2² (0} + 02)} Note the last operator is exactly the Laplacian. On the other hand, we have fundamental differential invariants x² + y² and u. To enforce symmetry, we replace the mamually defined set of extended derivatives {E*} in D-CIPHER by all nontrivial differential functions obtained from applying of searching for general functions of all variables (in this case, x, y, u) for the RHS expression, we restrict the search space to functions of fundamental differential invariants, i.e. x² + y² and u. Since D-CIPHER uses genetic programming to find a froe-form expression 9, we can simply replace the variable set in genetic programming by (1² + y², u to achieve this.
operator in O to one of the fundamental differential invariants. Also, instead
Table 6: Discovery results of D-CIPHER-based methods on the Darcy flow dataset. equation is 8(zuz yuy) (ugz + Ugy) +xx =0
The ground truth
Method D-CIPHER
Discovered equation
Dur + yuy 2.09xuy 2.09yu, 0.19u, 7.98x 82 2.51xy + 0.80 (xu + yuy,) 0.13 (Ux.r + Uyg) = 4.1
D-CIPHER-SI (ours)
Table shows the equations discovered by the D.CIPHER bascline and our method. Our method can find the correct functional form of the Darcy flow equation, while D-CIPHER with the original variables and derivative operators cannot. We comment that the benefit of symmetry is even greater here for D-CIPHER than for other SR methods like SINDy, because D-CIPHER requires the user to specify both the function coefficient a(x) and the function to be differentiated h(x, u) for an extended derivative. Such choices of functions can be largely arbitrary if no prior knowledge is available. On the other hand, our symmetry-based approach automatically selects this dictionary of differential functions.
C.6 Samples of Discovered Equations
In Table 7 we list some randomly selected equations discovered by different methods for the Boussinesq oquation Some mothods almost consistently discover correct/incorree equations (i.c., have success probabilitics close to or 0), so we only select one sample for each. For other methods with a large variance in the discovered equations, we display two samples: correct equation and incorrect one.
33

---

Published in Transactions on Machine Learning Research (03/2026)
The ground truth equation in the original variables is given in (6). The ground truth equation in the symmetry invariants is given by
(0.2) + ((00) + (4.0) + =
(50)
Table 7: Samples of discovered equations from the observed solution of the Boussinesq equation 6 For GP-based methods, include results from different numbers of itcrations (indicated by "N its"). For transformer-basee. methods, we include two samples for each method because of the large variance of discovered equations from different runs.
Method
Equation sample(s)
PySINDy PySINDy*
t1 -1 -1.01 0.79uuzz
Sparse
- 1.00 Myzza 0.99u - 0.98umaz
regression
S1 PySR (5 its) PySR (15 its) S1 (5 its)
(0.2)) = -1.00 - (4.0) - 1.00m7.00)72200
WUgz + 1.00uu + Ugzzr
Genetic
Wyr + Wet + 22 + 1.00mt...... =0
programming
1. 1.00moo.077(..00 + 1.007 + 1.007((4.0) + 11
(1) Uut - 1113uu 0.98u - 0..0 Un
E2E
(2) dy: -0.85uwz 0.75u 0.99u
Transformer
(1) M(,22) -.005 1.000)(4.0) 0.96
ST
(2) 70,2) -0.81710.0772.00 0.400 0.98m 0.90
Table 8 lists the equation samples discovered from the Darcy flow dataset. The ground truth equation in original variables is given in 77, and the ground truth equation in symmetry invariants is given by
8(2 Au - AR² 0,
(51)
where C2 = æug + y'ly, Au Ugw + Byy: and R² = x² + y² are among the rotational invariants used symbolic regression.
Table 8: Samples of discovered equations for the Darcy flow dataset.
Method
Equation sample
Genetic
PySR
0.47x2 0.38 0.09 +uyy) + 0.20 0
programming
SI
Ço 0.13Av
0.13e 1.01
E2E
Uzz 7.43 2² + 0.65u
Transformer
SI
Au = 2..66u + 0.85(2 + 0.29
Finally, Table 9 lists the equation samples discovered from the reaction-diffusion dataset. The ground truth
equation in original variables is given in (8) with di d2 0.1, and the ground truth equation in symmetry
invariants is given by
I. 0.1 + I + A(1 A)
E, = 0.1 + Eyy) A²
(52)
where Ip = uuy + Ven and E,- = -u + uw for any multiindex of t, and A u² u
C.7 Prediction Errors of Discovered Equations
In Table we report the prediction errors of the discovered equations on the three PDE systems. Specifically, for the Boussinesq equation and the reaction-diffusion system, we simulate the discovered PDE from an initial condition for a certain time period, c.g., E [0, 20) for the Boussinesq equation and t [0, 10] for the reaction-diffusion system. Then, we compare the numerical solution with the ground truth solution from the same initial condition at the end of the time period.
34

---

Published in Transactions on Machine Learning Research (03/2026)
Table 9: Samples of discovered equations for the reaction-diffusio system dataset. Each discovered result contains two equations, since this is an evolution system with two dependent variables u, v.
Method
Equation sample
Ut 0.96u 0.97u + 1.00 0.97mu + 1.00u2v + 0. 0.09ugy
PySINDY
Ut 0.96c 1.00u 0.97 1.00ue 0.96u² + 0.09 0.09k
Sparse regression
W₁ 0.21u 0.24u 1.00v
0.23uv² + 0.99u2v
PySINDy
= 0.21c - 1.01 - 0.24 - 0.99uv² 0.23u²
I. 0.101 + 0.101 vy + 0.96.1 0.96.4
SI
E, = 0.10E + 0.10E yy - 1.00.42
u, = 0.95u
0.96u + 1.00 0.96ur
1.00u² + 0.09u + 0.09uyv
SI-aligned
U2 0.95v
1.00 0.960
1.00uv
0.96u 0.09 + 0.09
iM = 0.92v
Genetic programming
PySR
6 -0.92
4 0.101
+ 0.101 + 1.00.4
E. 0.10
yz + 0.10 1.00.42
ut 0.89
E2E
U1 -0.9
Transformer
J,
ST
E 0.50 Aarctan(0.4552 - -5.1,, + +...) + +..))....
Prediction error for Boussinesq equation
Prediction error for reaction-diffusion system
10-3
SINDy-SI
SINDy-Sl-aligned
SINDy
SINDY-SI-raw
x 10
GP-SI
SINDY
GP
GP-SI
Transformer-Si
0.015
Transformer
RMSE
RMSE
0.010
10-2
0.005
4 10-2
0.000
30-1
0.0
2.5
5.0
7.5
10.0 Time
12.5
15.0
17.5
20.0
0
10
Figure 8: Prediction error over time using the discovered equations.
In addition to the prediction error the end of the simulation time, Figure shows the errors at each simulation timestep. We do not include methods whose error curves grow too fast due to the incorrectly identified equations. The results in Figure X are consistent with those in Table 1 Generally, the discovered equations with smaller prediction errors at the end of the sinulation time also have lower prediction errors throughout the entire time interval.
For Darcy flow 7 since describes the steady state system and does not involve time derivatives, we do simulate the discovered PDEs. Instead, evaluate each discovered PDE F(x, u(*) = 0 on the test dataset. (()) E 2. and report the residual the prediction error. Tn addition to the average error
35

---

Published in Transactions on Machine Learning Research (03/2026)
over all the spatial grid points reported in Table 1
visualize the error heatmaps the grid in Figure g
It can be observed that the discovered equations with symmetry invariants have lower errors across the entire grid.
Discovered equations with regular variables
Discovered equations with invariants
0.005
0.004
0.004
0.003
0.002
-00002
0.001
0.001
0.000
0.000
Figure 9: Prediction error of discovered equations from genetic programming methods for Darcy flow. Left: genetic programming with regular variables. Right: genetic programming with symmetry invariants.
D Comparison with Other Symmetry-Based Methods
In this section, discuss the connections and differences between work and other closely related methods that also enforce symmetry in equation discovery (Otto et al. 2023 Gurevich et al. 2024 Yang et al. 2024
D.1 Comparison with Methods for ODE Symmetries
Otto et al. (2023); Yang al. (2024) focus on & special type of Lie symmetry: time-independent (TI) symmetry of ODEs. In particular, the symmetry only transforms the phase variables of ODEs. In comparison, our method can handle equations with partial derivatives and symmetries acting on the independent variables, including time and spatial variables.
When restricted to the special case of ODEs linear combination forms and with TI symmetry, our method becomes equivalent EquivSINDy (Yang et al. 2024). More specifically, we can still follow the procedure in Section 4.2 to construct the invariants w.r.t the specified TI symmetry, and then apply Theorem 4.4 convert the symmel.ry constraint into linear constraints the SINDy parameters. This leads to the equivariant basis for the SINDy parameters in EquivSINDy.
Otto et al. 2023 Yang et al. h024 also introduce symmetry regularization, which is useful when computing the exact symmetry constraint challenging (for example, when symmetry is learned by a neural network instead of presented in closed form). The symmetry regularization term is based on the infinitesimal criterion of Lie point symmetrics, which applics Thus, the idea of symmetry regularization can be readily generalized to systems considered in our paper. paper primarily focuses enforcing hard symmetry constraints, choose not to investigate the effect of PDE symmetry regularization in great detail.
only ODEs but also PDEs and more complex symmetries.
However,
D.2 Comparison with SPIDER: Traaslaannnti Symmetry in 3D Fluid Systems
SPIDER (Gurevich et al. 2024 incorporates the symmetry rotations and translations in 3D space into sparse regression with weak formulation. With these specific designs, including symmetry, weak formulation,
36

---

Published in Transactions on Machine Learning Research (03/2026)
and sparse regression, their method is shown to successfully recover several canonical equations in fluid dynamics. Notably, their method for enforcing symmetry is similar to ours: they manually construct set of invariant scalars and a set of equivariant vectors and use them as feature sets for sparse regression.
Despite this similarity, we have established a more general framework for different types of symmetries. Examples in this paper include not only spatial rotation symmetry (Darcy flow 77), but also scaling symmetry
(Boussinesq equation (6) and phase-space symmetry (reaction-diffusinn { For each of these different
symmetries, can follow the standard procedure described in Section 4 to compute their invariants. Also, in terms of the base method, SPIDER primarily focuses the weak formulation of sparse regression. On the other hand, our method integrates with sparse regression (SINDy & WSINDy), genetic programming (PySR), and symbolic transformers.
To provide a direct experimental comparison to SPIDER (uurevich et al. 2024 we test our method (applied
SINDy) on the same channel flow dataset specified in thoir Table [Appendix A). The channel flow data
is retrieved from the Johns Hopkins Turbulence Database (Li al., 2008
Graham et al. 2016). The ground
truth equation is given by u, -u V)u Vp vAu, where V 5 10
The dataset contains a 3D velocity field and a scalar pressure field over 4D spatiotemporal grid (z, V. t). Our method requires set of scalar invariants of the given symmetry group. In this channel flow dataset, the assumed symmetry is rotations and spatiotemporal translations. To align with the assumption of STNDy, we specify the scalar invariant u u, as the LHS, and the following up-to-second-ores scalar invariants as the RHS variables:
(/ul2.p. V-n,ww n.w,
u.w. I-u V)u], Tr((²), Tr(s", /pl u V.V An, Ap, Vp-Au, w. w-Au, (An?), (53)
where S Vu+ (Vu)7] is the strain rate tensor, and w = V x u is the vorticity. We have excluded the invariants that involve mixed derivatives (such as Uay).
Table 10: Equation discovery results of SPIDER Gurevich et al. 2024 and our method on the channel flow dataset from the Johns Hopkins Turbulence Database.
Method
Discovered equation
SPIDER (vector)
1.000u; + 1.000(u V)u + 1.00VV - 0. 0.0000500Au SPIDER (scalar) u. 00 (u Vju 1.000u Vp + 0.0000473u Au 0. 0.497||u||=V
0.0000069|| Vull
SINDy-SI (ours)
u, 0.730V u 1.008u (u V)u 1.020u Vp 0.0000573u Auu
Table shows the discovered equations by our method (SINDy-ST) and SPIDER. The SPIDER results are referenced from their Table 1 and 2. It can be seen that SPIDER with the vector library (C1 in their paper) recovers the governing equation exactly. Our method and SPIDER using an extended scalar library (Co in their paper) discover mostly correct equations, with the pressure gradient, the convective term, and the Laplacian correctly identified. However, our method discovers one spurious term, the divergence V u, and SPIDER with the scalar library discovers two spurious terms (I|Vu||² and ||u|/²V- - u). While the divergence term should be zero in theory in the incompressible flow, the numerical simulation and derivative estimation may cause a small nonzero divergence, which then reflected in the equation discovery models.
We comment that, in this case, using the vector library in SPIDER is natural choice, since can easily obtain the equivariant vectors w.r.t rotations. SPIDER with the vector library also achieves the best accuracy. Still, we show that our framework, which requires an invariant scalar library, can also be applied to this scenario and identify a mostly correct equation with only one spurious term.
E Experiment Details
In this section, describe the experiment setups required to reproduce the experiments.
37

---

Published in Transactions on Machine Learning Research (03/2026)
E.1 Data generation
Boussinesq equation The equation is solved using a Fourier pseudospectral method for spatial derivatives and a fourtb-order Kutge-Kutta (RK4) scheme for time integration. The solution is computed on a periodic spatial domain [-10, 10] with N = 256 grid points. The equation is reformulated as a first-order system in time by introducing v = ue, and both u and v are evolved in time. Spatial derivatives are computed using the Fast Fourier Transform, and time derivatives of u up to the fourth order are derived analytically from the governing equation. At each time step, values of are recorded in the dataset for equation discovery. The simulation starts from initial condition of u(x) 0.5e - and u, 0 and proceeds up to final time T= 20 with time step of Al = 0.001. Starting from the solution at T = 20, we simulate for another T' = 20 with the same configuration to obtain a test dataset for evaluating prodiction errors of the discovered equations.
Darcy flow We use the data generation code from PDEBench (Takamoto et al. 2022 to generate the steady-state solution of Darcy flow over a unit square. The solution is obtained by numerically solving a temporal evolution equation
Ue (x t) V(((x)V f(x), (-0.0, (-5.5,0.5)
(54)
with a(x) = e- 4|x|| fx) = 1, a smooth random initial condition generated by the init init_multa multi 2DRand routine from PDEBench, and homogeneous Neumann boundary conditions (zero normal flux) on A. We integrate from t to & using an explicit two-stage scheme with an adaptive time step chosen from a diffusive CFL condition, with a CFL safety factor of 0.25.
We the data, generation code from PySINDy (de Silva et al. 2020 Kaptanoglu et al. 2022). We solve a periodic spatial domain of 10, 10] x [-10, 10 with a 128 128 Fourier spectral grid. The initial condition is given by u(x, y, 0) = tanh(r) cos(0 - r), v(x, y. 0) tanh(r) sin(O r), where V x/2 + y² and 0 arg(# + ig). The simulation uses RK45 in Fourier space and proceeds up to final time T = 10 with a time step At 0.05. We perturb the numerical solution by 0.05% noise and record the values of u, to the dataset for equation discovery. Starting from the solution at T = 10, we simulate for another T = 10 with the same configuration to obtain a test dataset for evaluating prediction errors of the discovered equations.
Reaction-diftusion
E.2 Sparse regression
Boussinesq equation For SINDy with original variables, we fix un as the LHS of the equation and include functions of up to 4th-order derivatives on the RHS. For PySINDy in Table [ the library contains monomials on (44 with degree in u no larger than and degree in any partial derivative terms ua no larger than 1. For example, u² u is included, but u², uz are not. For PySINDy", the library contains all monomials on (T(4) up to degree? 2. For example, uz and wu, are included. Note that the PySINDy* library does not contain all functions in the original PySINDy library, e.g., u² u²u, is not included because it has degree 3.
Our method, SI, uses the invariant set in Example A.6 for sparse regression. Specifically, To.22 = un/u is used as the LHS of the equation, and the rest of the invariants are included in the RHS. The function library contains all monomials of these RHS invariants up to degree 2. Also, since the invariants contain rational functions with u, on the denominator, we remove the data points with [ua] < 0.1 to avoid numerical issues.
We also conduct an additional experiment investigate the impact of the threshold value for |uz/- Tn Table 11 enumerate different threshold values from (0.0001, 0.001, 0.01, 0.1, 0.1, 0.3), and report the resulting filtered dataset sizes (and their proportions compared to the unfiltered dataset), and the success probability (SP) and the prediction error (PE) metrics as in Table [
First of all, we notice that when the threshold value is small (c= 0.0001), i.c. effectively no filtering, the success probability for SINDy using invariant functions dramatically decreases. This exactly shows the
httpp///dhhbuaammmpppeeee/////ep http:///wtlllllaooo////////////g 10 PDEFIND P
/data gen/data gen NLE/Reaction Diffusionn
38

---

Published in Transactions on Machine Learning Research (03/2026)
necessity of applying this numerical filter, ut values close to zero would cause the invariant features to have large magnitudes and make the SINDy optimization unstable.
Then, as we increase we observe that our method can achieve 100% success probability for e E {0.001, 0.01, 0.1}. showing its robustness to different choices of the threshold to some extent. When we further increase c, the filtered dataset becomes much smaller, and the success probability decreases. However, even with c c= 0.3 and only 99 data points, our method is still able recover the correct equation with more than 50% probability.
Table 11: SINDy with invariant functions on the Boussinesq equation when removing data points with |zl <c c for different threshold values c. In the second row, we report the number of samples in the filtered datasets and their proportions compared to the original dataset. The success probability (SP) and the prediction error (PF) are computed from 100 runs with different random seeds, in the same way Table 1 The prediction error is reported as median [25% quantile, 75% quantile).
Threshold
Dataset size
SP
PE
0.0001
99,756 (97.4%)
0.36
NaN
00.129, NaN
0.001
97,956 (95.7%)
1.00
0.103
[0.099, 0.118
0.01
85,591 (83.6%)
1.00
0.098
[0.098, 0.099
0.1
26,231 (25.6%)
1.00
0.098
[0.098, 0.098
0.2
1,318 (1.3%)
0.91
0.098
(0.097, 0.108
0.3
99 (0.1%)
0.52
0.100
[0.098, NaN]
For all methods, we flatten the data on the spatiotemporal grid and randomly sample 2% of the data for each run. The data filtering process in Sl-raw is performed after subsampling. The threshold value for sequential thresholding is set to 0.25, and the coefficient for L2 regularization is set to 0.05.
Darcy flow Sparse regression-based methods are not directly applicable to Darcy flow 7 because there exist terms such as +=(("++)²) While it still possible to include all necessary terms the function library so that the equation can be written in the linear combination form { the knowledge of these complicated terms is nontrivial and should not be assumed available before running the equation discovery algorithm.
Reaction-Diffusion For SINDY with original variables, We fix u, and is the LHS of the equation and include functions of up to 2nd-order spatial derivatives on the RHS. In PySINDy, the library contains monomials of up to degree 3 and all spatial derivatives to order 2. In PySINDy", the library contains all monomials of u, v and their up to second-order spatial derivatives up to degree 3.
Our method uses the invariant set (1, y, u² + U'u U(u uph. where u = uu, and a is a multiindex of t, xx y. We will denote 44 = u u, and E, u u,, We use 4, and E, as the LHS of the equation, and the rest of the invariants are included in the RHS. The function library contains all monomials of these RHS invariants up to degree 2.
We randomly sample 10% of the data for each run. The threshold value for sequential thresholding is set to 0.05. The coefficient for L2 regularization is set to for SINDy with original variables and 0.1 for our method with symmetry invariants.
For the experiments with different levels of noise (Section 5.4 we use weak SINDy as the base algorithm. We use the implementation of weak SINDy from the PySINDy package (Kaptanoglu et al. 2022). The function library is the same SINDy described above. To enforce symmetry, instead of directly using the symmetry invariants, we derive set of linear constraints on the sparse regression parameters to adapt to weak SINDy. This procedure is further described in Appendix B.3
E.3 Genetic Programming
In all experiments, determine an equation matches the ground truth we first expand the prediction into of monomial terms. We then eliminate all terms whose relative coefficient is below 0.01. For each
39

---

Published in Transactions on Machine Learning Research (03/2026)
term in the filtered expression, we see if it matches any term in the ground truth expression. This is done by randomly sampling 100 points from the standard normal distribution and evaluating both the prediction and candidate ground truth term on the generated points. Note that we drop the coefficients before evaluation. If all evaluations of the predicted term have a relative error of less than 5% from those of the ground truth, the terms are said to match. If there is a perfect matching between the terms in the ground truth and prediction, the prediction is listed as correct.
Rather than directly returning single equation, PySR finally produces hall-of-fame that consists of multiple candidate solutions with varying complexities. To finally pick single prediction, we use selection strategy equivalent to the "best" option from PySR.
Boussinesq equation For the Boussinesq equation ( we first randomly subsample 10000 datapoints. We configure PySR to use the addition and multiplication operators, to have 127 populations of size 27, and to have the default fraction-replaced coefficient of 0.00036.
When running with ordinary variables, we sequentially try fixing the LHS to cach variable in (x. u(() and allow the RHS to be a function of all remaining variables. Similarly, runs using invariants sequentially fix the LHS from the set given by Example A.6 and the RHS as function of all other invariants.
For each iteration count of 5, 10, and 15, we run the algorithm using invariant or ordinary variables and report the number of correct predictions out of 100 trials.
Darcy flow In the Darey experiment we eliminate all points that are within 3 pixels from the border and then randomly subsample 10000 datapoints. We configure PySR to use the addition, multiplication, and exponential operators; to have 127 populations of size 64; and to have a fraction-replaced coefficient of 0.1. We further constrain it to disallow nested exponentials (e.g. exp(exp(x) + +4)
We try all possible ordinary variables in (x, () for the LHS and the RHS is then a function the unused variables. Likewise when using invariants, we fix the LHS to each possible invariant specified in Example A.5 and set the RHS as a function of the remaining invariants.
For each iteration count of 50, 100, and 200, we run the algorithm using invariant or ordinary variables and report the number of correct predictions out of 100 trials.
Reaction-Diffusion For the Reaction Diffusion equation we remove all points that are within 3 pixels from the border or have timestamp greater than or equal to 40, and then randomly subsample 10000 datapoints. We configure PySR to use the addition and multiplication operators, have 127 populations of size 64, and to have fraction-replaced coefficient of 0.5.
In the ordinary variable case, we fix the LHS either uu eu and allow the RHS to be a function of all other variables in (x, (2) When using invariants, the LHS is fixed to be eit T, or F, and the RHS is then a function of all remaining invariants.
For each iteration count of 100, 200, and 400, we run the algorithm using regular and ordinary variables and report the number of correct predictions out of 100 trials.
E.4 Symbolic Transformer
We use the pretrained symbolic transformer model provided in the official codebase from Kamienny et al. (2022). The transformer- baseed symbolic regressor is initialized with 200 maximal input points and 100 expression trees to refine. The variable sets used in the symbolic transformer are the same as those described in the genetic programming experiments, except for the Boussinesq equation, where we remove all mixed derivative terms in both the original variable set and the symmetry invariant set. We find that the symbolic transformer can sometimes discover the correct equation under this further simplified setup, but fails when using the larger variable sets.
40

---

Published in Transactions on Machine Learning Research (03/2026)
We also fix the LHS of the function and use the remaining variables as RHS features. For the Boussinesq equation, the LHS is fixed to use for original variables and ((02)) for symmetry invariants. For the Darcy flow, the LHS is fixed to Uaz for original variables and Au for symmetry invariants. For the reaction-diffusion system, the LHS is fixed to ue, Ut for original variables and 4, Eg for symmetry invariants.
E.5 Hypothesis Spaces of Equation Discovery Algorithms
Table 12 and Table 13 describe the hypothesis spaces of different equation discovery algorithms when applied to the Boussinesq equation and the reaction-diffusin system.
Table 12: Hypothesis spaces of different equation discovery algorithms for the Boussinesq equation.
Method
Hypothesis space
PySINDy
Ver = ) (0) = {ab E Mono E {1, My, Ueara
Sparse
PySINDy*
Net P(u()) E =2 (U(4)
Regression
PySINDy"
Ue = ) (0') {u [u u uy us G 2 0. 2. 4 S Ei sgn(c;) S 1}
SI
[(.22) = P(n) E Polyca((V) (((o,))))) (Too.)))
Genetic PySR
z² = f(z +) for = (x, u (1) and some j
Programming
SI
(a0.80) = ((n) for n (M).
3 I 4] and some (an, 3o)
Table 13: Hypothesis spaces of different equation discovery algorithms for 2D reaction-difllsion. u(*) E U
denotes the collection of all up to nth order spatial derivatives. [a1, 02] is the multiindex for spatial variables. x (z, y. t). A = + +
Method
Hypothesis space
PySINDY PySINDy* PySINDy**
u; wo(u (00 Mu |a| 2 u, P(u()) (2) Poly <3 (U(2))
Sparse Regression
() ()), [00} = {I (w.) E < 3, 3, 0021 (g.) S 1
ST
14. E, P
(A,x, Eni |a| <
Sl-aligned
u, (w (2), Wik Quh33 8r for some precomputed Q
Genetic
PySR
u, f(x, (2))
Programming
S1
4, Eg] T f f(A,x,AA44 Las Eai a s 2)
41


Reviews:
[{"review_id": "odYnyek8BQ", "reviewer": "Reviewer_hkwD", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This article proposes a regression-based equation learning framework that relies on differential invariants of symmetry transformations as the fundamental building blocks, rather than raw variables. The goal of doing so is to ensure that the discovered equation satisfies existing physical symmetries.\n\nThis article is clear, well-written, and technically correct.  Most of the comments below are designed to improve the impat ofthe paper.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The authors propose their methodology and describe it in great detail. The efficacy of their proposed method is demonstrated on a handful of canonical examples. They do present some theory in the appendix, but this one lemma  is focused on the differential invariant symmetry transformations as the fundamental building blocks, rather than raw variables. The second does provide some comfort regarding convergence. However, convergence is different when using the strong vs. weak forms of SINDy, so it's not clear how translatable the theory is.\n\nIn support of the theory, then demonstrate successful learning of the equations for several PDE examples (even in the presence of noise).", "audience_yn": "Yes", "audience_explanation": "One of the major challenges of using equation error-based equation learning is that it is far too easy for non-physical terms to appear, and then the discovered model is completely wrong.  This article addresses that lack head on.", "requested_changes": "* Figure 3 SINDy -> WSINDy\n* Are the authors using WSINDy or the weak form implementation in pysindy. As WSINDy autotunes its internal parameters, it can be directly applied to data. The weak form implemented in pysindy is not the same as WSINDy and does not autotune parameters (like the test function radius).  There are several examples of Pysindy failing where WSINDy succeeded.\n* Given the strong similarity to the SPIDER paper by Gurevich et al., 2024 in JFM, the authors must compare with the results of that code. Of course, either the SPIDER code or the authors' code must be altered to be comparable.  The goals and overall approach of the two methods are just so similar that the authors need to include a serious discussion of the similarities and differences.\n* While there are several examples of implementing the proposed framework to enforce symmetries of differential equations. I would like to see a discussion on how easy it is to adapt an algorithm with the symmetry constraints described. They are described in technical detail, but a qualitative review would be helpful. Perhaps some of the material in Appendix C.3 can be summarized in 5.3.\n* In the paper, the authors state, “For more general symmetries and physical systems, enforcing symmetry often requires additional assumptions on the form of equations, such as the linear combination form in sparse regression (Otto et al., 2023; Yang et al., 2024).” One of my questions is how the proposed work compares to this method when an ODE can be written as a linear combination. Perhaps in general, it would be good to make a comparison to the other methods the authors reference that are able to enforce symmetry.", "broader_impact_concerns": "The broader impact is fine. I have no comments.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "kXjqZtL54B", "reviewer": "Reviewer_YyRe", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present a novel approach for symbolic discovery of differential equations by using symmetry invariant groups as the atomic entities in equation discovery. The approach is sophisticated but well-explained and involves learning a number of low-order invariants then composing those to construct sets of arbitrary orders. The paper is extremely well-written. The authors then show the behavior of their method on three PDE examples. The new method performs favorably to other methods such as SINDy and genetic programming.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The authors present results on three PDEs, and even explore the impact of noise. The results are convincing. That said, there are some issues (see below).", "audience_yn": "Yes", "audience_explanation": "This work is very timely in the field of scientific machine learning. Discovering equations from data in a fashion that respects underlying group structure will be critical for both trustworthy surrogate modeling and for learning in data-poor regimes.", "requested_changes": "I would recommend the following changes. These are crucial both to strengthening the work and for acceptance.\n\n1. For the PDEs that are presented, the authors should provide a complete description in the main body of the paper. It should include the domain in which the PDE is solved, the dimension, and the and the boundary conditions. Without this information, it's almost impossible to replicate results. That said, it looks like all the PDEs are in two spatial dimensions at most.\n\n2. I'd have liked to see both 3D and higher dimensional PDE examples. It's in these scenarios that discovery of symmetries will be particularly important, and these are the scenarios with greatest practical applicability.\n\n3. It is important to see nonlinear examples as well.\n\n4. Could the authors show with one example how the method performs when third or higher order invariants are involved? Perhaps an ablation over the order of invariants, if that's at all possible? Given enough solution smoothness, that should be doable in principle, though there may be practical difficulties.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "vw0yrKv7hr", "reviewer": "Reviewer_FLuz", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a general framework for enforcing symmetry constraints in symbolic regression (SR) methods used for discovering partial differential equations (PDEs). The key idea is to replace the standard set of variables and derivatives with differential invariants of a given symmetry group. This ensures that any equation discovered by the SR algorithm inherently respects the specified symmetry, thereby reducing the search space, improving parsimony, and preventing unphysical solutions.\n\nStrenghts:\n-The method is derived from the well-established theory of differential invariants and can be integrated with various existing SR backbones (sparse regression, genetic programming, transformers).\n-The paper provides comprehensive experiments on canonical PDE systems (Boussinesq, Darcy flow, reaction-diffusion), demonstrating improvements in accuracy, efficiency, and robustness—especially in challenging scenarios with noisy data and imperfect symmetry.\n-The proposed relaxation of strict symmetry constraints (Section 4.4) is a pragmatic solution for handling systems with approximate or broken symmetries, showing better performance than an unconstrained baseline.\n\nWeaknesses:\n-As already noted by the authors, the method assumes the symmetry group is known a priori. The framework cannot discover symmetry itself and may produce incorrect results if the symmetry is misspecified.\n-While the literature review is excellent within the specific scope of symbolic regression, its impact could be enhanced by briefly situating the work within the broader, compelling challenge of learning or discovering governing equations from data.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The claims made in the submission are supported by accurate, convincing, and clear evidence, as demonstrated through a multi-faceted validation strategy:\n\n-The core claim, that symmetric PDEs can be expressed via differential invariants, is anchored in established mathematical theory. The paper correctly leverages Theorem 4.2 as its foundation and provides derivations for constructing invariants.\n-The method is tested on three distinct, canonical PDE systems (Boussinesq, Darcy flow, reaction-diffusion) covering different challenges (high-order derivatives, multiple variables, steady-state vs. evolution).\n-Performance is compared fairly against standard implementations of three major SR algorithm classes (SINDy, PySR, E2E transformer) both with and without the symmetry invariant (SI) framework.\n-Success is measured using clear, relevant metrics: Success Probability (SP) for recovery of the exact symbolic form, and Prediction Error (PE) for the predictive fidelity of the discovered model.\n-Experiments systematically vary key conditions to test robustness: e.g. Noisy Data or Imperfect Symmetry\n-The integration of SI into different SR backbones is explained with explicit algorithms.\n- Implementation details, data generation procedures, and hyperparameters are thoroughly documented in the Appendices.\n-The paper acknowledges limitations openly.", "audience_yn": "Yes", "audience_explanation": "Yes. The paper’s integration of symmetry theory with scalable symbolic regression directly addresses a core interest in interpretable PDE discovery, a topic relevant to researchers in machine learning, scientific computing, and computational physics within the TMLR audience.", "requested_changes": "The paper is a well-executed and significant contribution that merits publication. My recommendation for acceptance is not contingent on any changes. However, the following minor suggestion would strengthen the work by providing a more comprehensive context for the broader TMLR audience.\n\nMinor Suggestion to Strengthen the Work:\nThe current related work section focuses primarily on symbolic regression methods. To better position this contribution within the extensive field of PDE discovery and learning from data, please consider adding a brief paragraph or subsection that acknowledges other major paradigms. This could include, but is not limited to, neural operator learning (e.g., DeepONet), Neural Networks and Gaussian process-based PDE learning, and methods for learning effective coarse-grained dynamics from fine-scale simulations. This addition would help readers understand how this symbolic, symmetry-informed approach complements and differs from non-symbolic or surrogate modeling techniques in the long-standing pursuit of discovering governing equations from data.\n\nThe current related work section focuses primarily on symbolic regression methods. To better position this contribution within the extensive field of PDE discovery and learning from data, I respectfully suggest the authors consider adding a brief paragraph or subsection that acknowledges other major paradigms. This could include neural operator learning, Gaussian process-based methods, and techniques for learning effective coarse-grained dynamics or closures from fine-scale simulations (a long-standing goal in numerical analysis and scientific machine learning).\n\nThe authors may of course select their own preferred set of references to represent these directions. For their consideration, the following works exemplify the breadth of related approaches:\n\n\nLee et al. (2020). Coarse-scale PDEs from fine-scale observations via machine learning. Chaos.\n\nVlachas et al. (2020). Backpropagation algorithms and reservoir computing for spatiotemporal dynamics. Neural Networks.\n\nLi et al. (2020). Hamiltonian neural networks. Proceedings of the 2nd Symposium on Advances in Approximate Bayesian Inference.\n\nChen et al. (2021). Solving and learning nonlinear PDEs with Gaussian processes. Journal of Computational Physics.\n\nLu et al. (2021). Learning nonlinear operators via DeepONet. Nature Machine Intelligence.\n\nGalaris et al. (2022). Numerical bifurcation analysis of PDEs from lattice Boltzmann model simulations. Journal of Scientific Computing.\n\nVlachas et al. (2022). Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelligence.\n\nFloryan & Graham (2022). Data-driven discovery of intrinsic dynamics. Nature Machine Intelligence.\n\nLee et al. (2023). Learning chemotactic PDEs/closures from agent-based data. Journal of Mathematical Biology.\n\nDietrich et al. (2023). Learning effective stochastic differential equations. Chaos.\n\nFabiani et al. (2024). Task-oriented machine learning surrogates for tipping points of agent-based models. Nature Communications.\n\nA concise discussion situating the current symmetry-informed symbolic approach alongside such non-symbolic or surrogate modeling techniques would be valuable for readers.", "broader_impact_concerns": "No broader impact concerns. The authors have included a clear and appropriate broader impact statement, acknowledging both the potential scientific benefits and the essential need for domain-expert validation to mitigate risks from erroneous equations. The work presents no identifiable or specific ethical risks.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)