Discovering Symbolic Differential Equations with Symmetry Invariants

## Abstract

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 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 fluid and reaction-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 a wide range of phenomena, from fluid dynamics and electromagnetic fields to chemical reactions and biological processes, as they succinctly capture the underlying principles governing the behavior of complex systems. The discovery of governing equations in symbolic forms from observational 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, we aim to discover symbolic partial differential equations (PDEs) in the form

$$ $F(\mathbf{x},u^{(n)})=0,$ (1) $$

where $\mathbf{x}$ denotes the independent variables, $u^{(n)}$ consists of the dependent variable $u$ and all of its up-to-$n$th order partial derivatives.

While it has long been an exclusive task for human experts to identify governing equations, symbolic regression (SR) has emerged as an increasingly popular approach to automate the discovery.(^1^11While some literature uses symbolic regression specifically for GP-based methods, we use the term interchangeably with equation discovery to refer to all algorithms for learning symbolic equations.) SR constructs expressions from a predefined set of atomic entities, such as variables, constants, and mathematical operators, and fits the expressions to data by numerical optimization. Common methods include sparse regression (Brunton et al., 2016; Champion et al., 2019), genetic programming (Cranmer et al., 2019, 2020; Cranmer, 2023), neural networks (Kamienny et al., 2022), etc.

However, symbolic regression 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 expedites the search process, and it also rules out physically invalid or unnecessarily complex equations.

Figure: Figure 1: 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 a circular domain and demonstrate that using symmetry invariants guarantees a symmetric output.
Refer to caption: extracted/6449288/figures/overview.png

Among the various physical constraints,
symmetry plays a 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 as 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) deals 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 symmetry invariants to enforce the inductive bias of symmetry with minimal restrictions in the types of equations and SR algorithms.
Specifically, we leverage the fact that a 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 in symbolic regression, as depicted in [Figure 1](https://arxiv.org/html/2505.12083v1#S1.F1). These 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-based 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_{t}\coloneqq\partial u/\partial t$, and $u_{xx}\coloneqq\partial^{2}u/\partial x^{2}$. 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, respectively, e.g. $W,\mathbf{w},w$. These conventions may admit exceptions for clarity or context. See [Table 2](https://arxiv.org/html/2505.12083v1#A2.T2) for a full description of notations.

## 2 Background

### 2.1 PDE Symmetry

This section introduces the basic concepts about 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 Differential Equations.

We consider PDEs in the form $F(\mathbf{x},u^{(n)})=0$, as given in ([1](https://arxiv.org/html/2505.12083v1#S1.E1)). We restrict ourselves to a single equation and a single dependent variable, though generalization is possible.
We use $\mathbf{x}\in X\subset\mathbb{R}^{p}$ to denote all independent variables. For example, $\mathbf{x}=(t,x)$ for a system evolving in 1D space. Note that the bold $\mathbf{x}$ refers to the collection of all independent variables while the regular $x$ denotes the spatial variable. Then, $u=u(\mathbf{x})\in U\subset\mathbb{R}$ is the dependent variable; $u^{(n)}=(u,u_{x},...)$ denotes all up to $n$th-order partial derivatives of $u$; $(\mathbf{x},u^{(n)})\in M^{(n)}\subset X\times U^{(n)}$, where $M^{(n)}$ is the $n$th order jet space of the total space $X\times U$. $M^{(n)}$ and $u^{(n)}$ are also known as the $n$th-order prolongation of $X\times U$ and $u$, respectively.

#### Symmetry of a PDE.

A point symmetry $g$ is a local diffeomorphism on the total space $E=X\times U$:

$$ $g\cdot(\mathbf{x},u)=(\tilde{\mathbf{x}}(\mathbf{x},u),\tilde{u}(\mathbf{x},u)),$ (2) $$

where $\tilde{\mathbf{x}}$ and $\tilde{u}$ are functions of $E$.
The action of $g$ on the function $u(\mathbf{x})$ is induced from ([2](https://arxiv.org/html/2505.12083v1#S2.E2)) by applying it to the graph of $u:X\to U$. Specifically, denote the domain of $u$ as $\Omega\subset X$ and its graph as $\Gamma_{u}=\{(\mathbf{x},u(\mathbf{x})):\mathbf{x}\in\Omega\}$. The group element $g$ transforms the graph $\Gamma_{u}$ as
$\tilde{\Gamma}_{u}\coloneqq g\cdot\Gamma_{u}=\{(\tilde{\mathbf{x}},\tilde{u})=
g\cdot(\mathbf{x},u):(\mathbf{x},u)\in\Gamma_{u}\}$.

Since $g$ transforms both independent and dependent variables, $\tilde{\Gamma}_{u}$ does not necessarily correspond to the graph of any single-valued function. Nevertheless, by suitably shrinking the domain $\Omega_{X}$, we can ensure that the transformations close to the identity transform $\Gamma_{u}$ to the graph of another function. This function with the transformed graph $\tilde{\Gamma}_{u}$ is then defined to be the transformed function of the original solution $u$, i.e. $g\cdot u=\tilde{u}$ s.t. $\Gamma_{\tilde{u}}=\tilde{\Gamma}_{u}$. The symmetry of the PDE $\eqref{eq:pde-def}$ is then defined:

###### Definition 2.1 .

A symmetry group of $F(\mathbf{x},u^{(n)})=0$ is a local group of transformations $G$ acting on an open subset of the total space $X\times U$ such that, for any solution $u$ to $F=0$ and any $g\in G$, the function $\tilde{u}=(g\cdot u)(\mathbf{x})$ is also a solution of $F=0$ wherever it is defined.

###### Definition 2.1 .

#### Infinitesimal Generators.

Often, the symmetry group of a PDE is a continuous Lie group.
In practice, one needs to compute with infinitesimal generators of continuous symmetries, i.e., vector fields.
In more detail, we will write vector fields $\mathbf{v}:E\to TE$ on $E=X\times U$ as

$$ $\mathbf{v}=\xi^{j}(\mathbf{x},u)\frac{\partial}{\partial x^{j}}+\phi(\mathbf{x },u)\frac{\partial}{\partial u}.$ (3) $$

Any such vector field generates a one-parameter group of symmetries of the total space {$\exp(\epsilon\mathbf{v}):\epsilon\in\mathbb{R}\}$. The symmetries arising from the exponentiation of a vector field moves a point in the total space along the directions given by the vector field. We will specify symmetries by vector fields in the following sections. For instance, $\mathbf{v}=x\partial_{y}-y\partial_{x}$ represents the rotation in $(x,y)$-plane; $\mathbf{v}=\partial_{t}$ corresponds to time translation.

Since we deal with PDEs, we need to consider the prolonged group actions and infinitesimal actions on the $n$th-order jet space, induced from ([2](https://arxiv.org/html/2505.12083v1#S2.E2)) and ([3](https://arxiv.org/html/2505.12083v1#S2.E3)). These are denoted $g^{(n)}$ and $\mathbf{v}^{(n)}$, respectively. A more detailed discussion on prolongation of group actions is deferred to [Section B.2](https://arxiv.org/html/2505.12083v1#A2.SS2). To introduce our method, it suffices to note that the prolongation of the vector field ([3](https://arxiv.org/html/2505.12083v1#S2.E3)) can be described explicitly by $\xi^{j}$ and $\phi$ and their derivatives via the prolongation formula ([11](https://arxiv.org/html/2505.12083v1#A2.E11)).

### 2.2 Symbolic Regression Algorithms

Given the data $\{(x^{i},y^{i})\}\subset X\times Y$, the objective of symbolic regression is to find a symbolic expression for the function $y=f(x)$. Although this original formulation is for algebraic equations, it can be generalized to differential equations like ([1](https://arxiv.org/html/2505.12083v1#S1.E1)). To discover a PDE from the dataset of its observed solutions on a grid $\Omega$, i.e., $\{(\mathbf{x},u(\mathbf{x})):\mathbf{x}\in\Omega\}$, we estimate the partial derivative terms and add them to the dataset: $\{(\mathbf{x},u^{(n)}):\mathbf{x}\in\Omega\}$. One of the variables in the variable set $(\mathbf{x},u^{(n)})$ 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 or speculations about the underlying system.

We briefly review two classes of symbolic regression algorithms: sparse regression (SINDy) and genetic programming (GP). A more detailed discussion of related works is found in [Appendix A](https://arxiv.org/html/2505.12083v1#A1).

Sparse regression (Brunton et al., 2016) is specifically designed for discovering differential equations. It assumes the LHS $\ell$ of the equation is a fixed term, e.g. $\ell=u_{t}$, and the RHS of the equation can be written as a linear combination of $m$ predefined functions $\theta^{j}$ with trainable coefficients $\mathbf{w}\in\mathbb{R}^{m}$, i.e.,

$$ $\ell(\mathbf{x},u^{(n)})=w^{j}\theta^{j}(\mathbf{x},u^{(n)}),\ \theta^{j}:M^{( n)}\to\mathbb{R}.$ (4) $$

The equation is found by solving for $\mathbf{w}$ that minimizes the objective
$\|\mathbf{L}-\mathbf{R}\|_{2}^{2}+\lambda\|\mathbf{w}\|_{0},$
where $\mathbf{L}$ and $\mathbf{R}$ are obtained by evaluating $\ell$ and $w^{j}\theta^{j}$ on all data points and concatenating them into column vectors, and $\|\mathbf{w}\|_{0}$ regularizes the number of nonzero terms. This formulation can be easily extended to $q$ equations and dependent variables ($q>1$): $\ell^{i}(\mathbf{x},\mathbf{u}^{(n)})=W^{ij}\theta^{j}(\mathbf{x},\mathbf{u}^{
(n)})$, $W\in\mathbb{R}^{q\times m}$.

One problem with sparse regression is that its assumptions about the form of the equation are restrictive. Many equations cannot be expressed in the form of ([4](https://arxiv.org/html/2505.12083v1#S2.E4)), e.g. $y=\frac{1}{x+a}$ where $a$ could be any constant. Also, the success of sparse regression relies on the proper choice of the predefined function library $\{\theta^{j}\}$. If any term in the true equation were not included, sparse regression would fail to identify the correct equation.

Genetic programming
offers an alternative solution for equation discovery (Cranmer, 2023), which is capable of learning equations in more general forms. It represents each expression as a tree and instantiates a population of individual expressions. At each iteration, it randomly samples a subset of expressions and selects one of the expressions 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 a weaker 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. Genetic programming can be less efficient than sparse regression when the equation can be expressed in the form ([4](https://arxiv.org/html/2505.12083v1#S2.E4)) due to its larger search space. However, we will show that it is a promising alternative to discover PDEs of generic forms, and our approach further boosts its efficiency.

## 3 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 methods.

### 3.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 3.1 (Def 2.51, Olver ( 1993 ) ) .

Let $G$ be a local group of transformations acting on $X\times U$. Any $g\in G$ gives a prolonged group action $\mathrm{pr}^{(n)}g$ on the jet space $M^{(n)}\subset X\times U^{(n)}$. An $n$th order differential invariant of $G$ is a smooth function $\eta:M^{(n)}\to\mathbb{R}$, such that for all $g\in G$ and all $(\mathbf{x},u^{(n)})\in M^{(n)}$, $\eta(g^{(n)}\cdot(\mathbf{x},u^{(n)}))=\eta(\mathbf{x},u^{(n)})$ whenever $g^{(n)}\cdot(\mathbf{x},u^{(n)})$ 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$ is generated by a set of infinitesimal generators $\mathcal{B}=\{\mathbf{v}_{a}\}$, then a function $\eta$ is a differential invariant of $G$ iff $\mathbf{v}_{a}^{(n)}(\eta)=0$ for all $\mathbf{v}_{a}\in\mathcal{B}$.

The following theorem guarantees that any differential equation admitting a symmetry group can be expressed solely in terms of the group invariants.

###### Theorem 3.2 (Prop 2.56, Olver ( 1993 ) ) .

Let $G$ be a local group of transformations acting on $X\times U$. Let $\{\eta^{1}(\mathbf{x},u^{(n)}),...,\eta^{k}(\mathbf{x},u^{(n)})\}$ be a complete set of functionally independent $n$th-order differential invariants of $G$. An $n$th-order differential equation ([1](https://arxiv.org/html/2505.12083v1#S1.E1)) admits $G$ as a symmetry group if and only if it is equivalent to an equation of the form

$$ $\tilde{F}(\eta^{1},...,\eta^{k})=0.$ (5) $$

Consequently, symbolic regression with a complete set of invariants precisely searches within the space of all symmetric differential equations, while automatically excluding equations violating the specified symmetry.

Our strategy of using differential invariants applies broadly to various equation discovery algorithms. For instance, in sparse regression, we can construct the function library using invariants rather than raw variables 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 3.2](https://arxiv.org/html/2505.12083v1#S3.SS2)), and how to incorporate them into specific SR algorithms ([Section 3.3](https://arxiv.org/html/2505.12083v1#S3.SS3)).

###### Definition 3.1 (Def 2.51, Olver ( 1993 ) ) .

###### Theorem 3.2 (Prop 2.56, Olver ( 1993 ) ) .

### 3.2 Constructing a Complete Set of Invariants

Despite the simplicity of our strategy, we still need a concrete method for computing the invariants.
In this subsection, we provide a general guideline to construct a complete set of differential invariants up to a required order given the group action.

By definition of differential invariants, we look for functions $\eta(\mathbf{x},u^{(n)})$ satisfying $\mathbf{v}^{(n)}(\eta)=0$ given a prolonged vector field $\mathbf{v}^{(n)}$. This is a first-order linear PDE that can be solved by the method of characteristics. However, in practice, if $E=X\times U\simeq\mathbb{R}^{p}\times\mathbb{R}$, there are $\binom{p+n-1}{n}$ partial derivatives of the independent variable $u$ of order exactly $n$. Therefore, as $n$ grows, it quickly becomes impractical to solve directly for $n$th-order differential invariants. The
computation of higher-order differential invariants, if necessary, is made tractable
by the following result, where higher-order invariants are computed recursively from lower-order ones.

###### Proposition 3.3 .

Let $G$ be a local group of transformations acting on $X\times U\simeq\mathbb{R}^{p}\times\mathbb{R}$. Let $\eta^{1},\eta^{2},\cdots,\eta^{p}$ be any $p$ differential invariants of $G$ whose horizontal Jacobian $J=[D_{i}\eta^{j}]$ is non-degenerate on an open subset $\Omega\subset M^{(n)}$. If there are a maximal number of independent, strictly $n$th-order differential invariants $\zeta^{1},\cdots,\zeta^{q_{n}}$, $q_{n}=\binom{p+n-1}{n}$, then the following set contains a complete set of independent, strictly $(n+1)$th-order differential invariants defined on $\Omega$:

$$ $\frac{\mathrm{det}(D_{i}\tilde{\eta}^{j}_{(k,k^{\prime})})}{\mathrm{det}(D_{i} \eta^{j})},\ \forall k\in[p],k^{\prime}\in[q_{n}],$ (6) $$

where $i,j\in[p]$ are matrix indices, $D_{i}$ denotes the total derivative w.r.t $i$-th independent variable and $\tilde{\eta}^{j}_{(k,k^{\prime})}=[\eta^{1},...,\eta^{k-1},\zeta^{k^{\prime}},
\eta^{k+1},...,\eta^{p}]$.

In practice, we first solve for $\mathrm{pr}\ \mathbf{v}(\eta)=0$ to obtain a sufficient number of lower-order invariants as required in [Proposition 3.3](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem3), starting from which we can construct complete sets of invariants of arbitrary orders. Direct results from applying ([6](https://arxiv.org/html/2505.12083v1#S3.E6)) can be complicated, especially for higher-order invariants. Thus, we often combine them to get simpler invariants, which we later use as the variable set in equation discovery.
In [Section B.4](https://arxiv.org/html/2505.12083v1#A2.SS4), we provide two examples of different symmetry groups and their differential invariants. Those results will also be used in our experiments.

###### Proposition 3.3 .

### 3.3 Implementation in SR Algorithms

Our symmetry principle characterizes a subspace of all equations with a given symmetry. Generally, this subspace partially overlaps with the hypothesis spaces of symbolic regression algorithms, conceptually visualized in [Figure 2](https://arxiv.org/html/2505.12083v1#S3.F2). As in [Theorem 3.2](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem2), PDEs with symmetry can be expressed as implicit functions of all differential invariants. However, symbolic regression methods typically learn explicit functions mapping features to labels. Some algorithms, such as SINDy, impose even stronger constraints on equation forms. Therefore, adaptation is needed to implement our strategy of using differential invariants in specific symbolic regression algorithms, as detailed below.

Figure: Figure 2: Venn diagram of hypothesis spaces from base SR methods and our symmetry principle.
Refer to caption: x1.png

#### General explicit SR

We start with general SR methods that learn an explicit function $y=f(x)$ without additional assumptions about the form of $f$, e.g., genetic programming and symbolic transformer. When learning the equation in terms of differential invariants, we do not know which one of them should be used as the LHS of the equation, i.e., the label $y$ in symbolic regression. Thus, we fit an equation for each invariant as LHS and choose the equation with the lowest data error, as described in [Algorithm 1](https://arxiv.org/html/2505.12083v1#alg1). We use the relative error to select the best equation because the scales of LHS terms differ.

Figure: Algorithm 1 General explicit SR for differential equations with symmetry invariants

#### Sparse regression

SINDy assumes a linear equation form ([4](https://arxiv.org/html/2505.12083v1#S2.E4)). Generally, its function library differs from the set of differential invariants. Also, SINDy fixes a LHS term, while we do not single out an invariant as the LHS of the equation when constructing the set of invariants.

Assume we are provided the SINDy configuration, i.e., the LHS term $\ell$ and the function library $\{\theta^{j}\}$. To implement sparse regression with symmetry invariants, we assign an invariant $\eta^{k}$ that symbolically depends on $\ell$, i.e., $\partial\eta^{k}/\partial\ell\not=0$, as the LHS for the equation in terms of symmetry invariants. For example, if $\ell=u_{tt}$ and the set of invariants is given by ([43](https://arxiv.org/html/2505.12083v1#A2.E43)), we use $\eta_{(0,2)}=u_{tt}u_{x}^{(b-2)/(a-b)}$ as the LHS since it is the only invariant that involves $u_{tt}$. 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 $\eta^{k}=\tilde{w}^{j}\theta^{j}(\bm{\eta}^{-k})$. Similar to [Algorithm 1](https://arxiv.org/html/2505.12083v1#alg1), we can expand all $\eta$ variables to obtain the equation in original jet variables.

The above approach optimizes an unconstrained coefficient vector $\tilde{\mathbf{w}}$ for functions of symmetry invariants. Alternatively, we can use the original SINDy equation form ([4](https://arxiv.org/html/2505.12083v1#S2.E4)) and implement the symmetry constraint as a constraint on the coefficient $\mathbf{w}$, as demonstrated in the following theorem. Here, we generalize the setup to multiple dependent variables and equations.

###### Proposition 3.4 .

Let $\bm{\ell}(\mathbf{x},\mathbf{u}^{(n)})=W\bm{\theta}(\mathbf{x},\mathbf{u}^{(n)})$ be a system of $q$ differential equations admitting a symmetry group $G$, where $\mathbf{x}\in\mathbb{R}^{p}$, $\mathbf{u}\in\mathbb{R}^{q}$, $\bm{\theta}\in\mathbb{R}^{m}$. Assume that there exist some $n$th-order invariants of $G$, $\eta_{0}^{1:q}\text{ and }\eta^{1:K}$, s.t. (1) the system of differential equations can be expressed as $\bm{\eta}_{0}=W^{\prime}\bm{\theta}^{\prime}(\bm{\eta})$, where $\bm{\eta}_{0}=[\eta_{0}^{1:q}]$ and $\bm{\eta}=[\eta^{1:K}]$, and (2) $\eta_{0}^{i}=T^{ijk}\theta^{k}\ell^{j}$ and $(\theta^{\prime})^{i}=S^{ij}\theta^{j}$, for some library functions $\bm{\theta}^{\prime}(\bm{\eta})$ and some constant tensors $W^{\prime}$, $T$ and $S$. Then, the space of all possible $W$ is a linear subspace of $\mathbb{R}^{q\times m}$.

Intuitively, the conditions in [Proposition 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4) state that the equations can be expressed as a linear combination of invariant terms, similar to the form in ([4](https://arxiv.org/html/2505.12083v1#S2.E4)) w.r.t original jet variables. Also, every invariant term in $\bm{\eta}_{0}$ and $\bm{\theta}^{\prime}(\bm{\eta})$ is already encoded in the original library $\bm{\theta}$. In practice, we need to choose a suitable set of invariants according to the SINDy configuration to meet these conditions. For example, when $\bm{\theta}$ contains all monomials on $M^{(n)}$ up to some degree, we can choose a set of invariants where each invariant is a polynomial on $M^{(n)}$. The proof of [Proposition 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4) is deferred to [Appendix B](https://arxiv.org/html/2505.12083v1#A2), where we explicitly identify the linear subspace for $W$ entailed by the proposition.

[Proposition 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4) allows us to keep track of the original SINDy parameters $W$ during optimization. This enables straightforward integration of symmetry constraints to variants of SINDy, e.g. Weak SINDy (Messenger & Bortz, 2021a, b) for noisy data. For example, if the constrained subspace has a basis $Q\in\mathbb{R}^{r\times q\times m}$, where $r$ is the subspace dimension, we write $W^{jk}=Q^{ijk}\beta^{i}$. While we directly optimize $\beta$, 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 take the form $\bm{\eta}_{0}=W^{\prime}\bm{\theta}^{\prime}(\bm{\eta})$, it is challenging to formulate the objective of Weak SINDy w.r.t $W^{\prime}$.

More implementation details related to [Section 3.3](https://arxiv.org/html/2505.12083v1#S3.SS3) can be found in [Appendix C](https://arxiv.org/html/2505.12083v1#A3).

###### Proposition 3.4 .

### 3.4 Constraint Relaxation for Systems with Imperfect Symmetry

Our approach discovers PDEs assuming perfect symmetry. However, it is common in reality that a system exhibits imperfect symmetry due to external forces, boundary conditions, etc. (Wang et al., 2022). In these cases, the previously mentioned method would fail to identify any symmetry-breaking 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 $\bm{\ell}=W\bm{\theta}$. We adopt the technique from Residual Pathway Prior (RPP) (Finzi et al., 2021), which is originally developed for equivariant linear layers in neural nets. 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^{jk}=Q^{ijk}\beta^{i}$ and $B^{jk}=P^{ijk}\gamma^{i}$ and place a stronger regularization on $\gamma$ than on $\beta$. While the model still favors equations in the symmetry subspace spanned by $Q$, symmetry-breaking components in $P$ can appear if it fits the data well.

## 4 Experiments

### 4.1 Datasets and Their Symmetries

We consider the following PDE systems, which cover different challenges in PDE discovery, such as 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 [Section E.1](https://arxiv.org/html/2505.12083v1#A5.SS1).

Water Wave.
Consider the Boussinesq equation describing the unidirectional propagation of a solitary wave in shallow water (Newell, 1985):

$$ $u_{tt}+uu_{xx}+u_{x}^{2}+u_{xxxx}=0$ (7) $$

This equation has a scaling symmetry $\mathbf{v}_{1}=2t\partial_{t}+x\partial_{x}-2u\partial_{u}$ and the translation symmetries in space and time. As shown in [Section B.4](https://arxiv.org/html/2505.12083v1#A2.SS4), the differential invariants are given by
$\eta_{(\alpha,\beta)}=u_{x^{(\alpha)}t^{(\beta)}}u_{x}^{-(2+\alpha+2\beta)/3}$
where $\alpha$ and $\beta$ are the orders of partial derivatives in $x$ and $t$, respectively. To discover the $4$th-order equation, we compute all $\eta_{(\alpha,\beta)}$ for $0\leq\alpha+\beta\leq 4$, except for $\eta_{(1,0)}=1$.

Darcy Flow.
The following PDE describes the steady state of a 2D Darcy flow (Takamoto et al., 2022) with spatially varying viscosity $a(x,y)=e^{-4(x^{2}+y^{2})}$ and a constant force term $f(x)=1$:

$$ $-\nabla(e^{-4(x^{2}+y^{2})}\nabla u)=1$ (8) $$

This equation admits an $\mathrm{SO}(2)$ rotation symmetry $\mathbf{v}=y\partial_{x}-x\partial_{y}$. A detailed calculation of the differential invariants of this group can be found in [Example B.5](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem5). In our experiment, we use the following complete set of $2$nd-order invariants: $\{\frac{1}{2}(x^{2}+y^{2}),u,xu_{y}-yu_{x},xu_{x}+yu_{y},u_{xx}+u_{yy},u_{xx}^
{2}+2u_{xy}^{2}+u_{yy}^{2},x^{2}u_{xx}+y^{2}u_{yy}+2xyu_{xy}\}$.

Reaction-Diffusion.
We consider the following system of PDEs from Champion et al. (2019):

$$ $\displaystyle u_{t}$ $\displaystyle=d_{1}\nabla^{2}u+(1-u^{2}-v^{2})u+(u^{2}+v^{2})v$ (9) $\displaystyle v_{t}$ $\displaystyle=d_{2}\nabla^{2}v-(u^{2}+v^{2})u+(1-u^{2}-v^{2})v$ (10) $$

In the default setup, we use $d_{1}=d_{2}=0.1$. The system then exhibits rotational symmetry in the phase space: $\mathbf{v}=u\partial_{v}-v\partial_{u}$. The ordinary invariants are $\{t,x,y,u^{2}+v^{2}\}$. The higher-order invariants are $\{\mathbf{u}\cdot\mathbf{u}_{\mu},\mathbf{u}^{\perp}\cdot\mathbf{u}_{\mu}\}$, where $\mathbf{u}=(u,v)^{T}$ and $\mu$ is any multi-index of $t$, $x$ and $y$.

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: $d_{1}=0.1$, $d_{2}=0.1+\epsilon$. 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 ([10](https://arxiv.org/html/2505.12083v1#S4.E10))is modified by adding $-\epsilon v$ to the RHS of $u_{t}$ and $-\epsilon u$ to the RHS of $v_{t}$. This can reflect a weak parametric forcing on the system.

### 4.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. (2022)), 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., $(\mathbf{x},u^{(n)})$ ) and our method using symmetry invariants. Our method will be referenced as SI (Symmetry Invariants) in the results.

To evaluate an equation discovery algorithm, we run it 100 times with randomly sampled data subsets and randomly initialized models if applicable. We record its success probability (SP) of discovering the correct equation. Specifically, we expand the ground truth equation into $\sum_{i}c^{i}f^{i}(\mathbf{z})=0$, where $c^{i}$ are nonzero coefficients, $\mathbf{z}$ denotes the variables involved in the algorithm, i.e., original jet variables $(\mathbf{x},u^{(n)})$ for baselines and symmetry invariants for our method, and $f^{i}$ are functions of $\mathbf{z}$. Also, the discovered equation is expanded as $\sum_{i}\hat{c}^{i}\hat{f}^{i}(\mathbf{z})=0$, where $\hat{c}^{i}\not=0$. The discovered equation is considered correct if all the terms with nonzero coefficients match the ground truth, i.e., $\{f^{i}\}=\{\hat{f}^{i}\}$.
We also report the prediction error (PE), which measures how well the discovered equation fits the data. For evolution equations, we simulate each discovered equation from an initial condition and measure its difference from the ground truth solution at a specific timestep in terms of root mean square error (RMSE). Otherwise, we just report the RMSE of the discovered equation evaluated on all test data points.

### 4.3 Results on Clean Data with Perfect Symmetry

**Table 1: Equation discovery results on clean data. $\mathcal{C}$, standing for complexity, refers to the effective parameter space dimension in sparse regression and the number of variables in GP/Transformer. SP and PE stands for success probability and prediction error, as explained in [Section 4.2](https://arxiv.org/html/2505.12083v1#S4.SS2). The entries "-" suggest that the method does not apply to the specific PDE system, or the result is not meaningful.**
| Method | Boussinesq ([7](https://arxiv.org/html/2505.12083v1#S4.E7)) | Darcy flow ([8](https://arxiv.org/html/2505.12083v1#S4.E8)) | Reaction-diffusion ([10](https://arxiv.org/html/2505.12083v1#S4.E10)) |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\mathcal{C}$ $\downarrow$ | SP $\uparrow$ | PE $\downarrow$ | $\mathcal{C}$ $\downarrow$ | SP $\uparrow$ | PE $\downarrow$ | $\mathcal{C}$ $\downarrow$ | SP $\uparrow$ | PE $\downarrow$ |  |  |
| Sparse Regression | PySINDy | 15 | 0.00 | 0.373 | - | - | - | 38 | 0.53 | 0.021 |
| SI | 13 | 1.00 | 0.098 | - | - | - | 28 | 0.54 | 0.008 |  |
| Genetic Programming | PySR | 17 | 0.90 | 0.098 | 8 | 0.00 | 0.187 | 17 | 0.00 | - |
| SI | 14 | 1.00 | 0.098 | 7 | 0.79 | 0.051 | 16 | 0.81 | 0.023 |  |
| Transformer | E2E | 10 | 0.53 | 0.132 | 8 | 0.00 | - | 17 | 0.00 | - |
| SI | 7 | 0.85 | 0.104 | 7 | 0.00 | - | 16 | 0.00 | - |  |

[Table 1](https://arxiv.org/html/2505.12083v1#S4.T1) summarizes the performance of all methods on the three PDE systems. For prediction errors (PE), we report the median, instead of the average, of 100 runs for each algorithm, because some incorrectly discovered equations yield tremendous prediction errors. Comparisons are made within each class of methods. Generally, using symmetry invariants reduces the complexity of equation discovery and improves the chance of finding the correct equations compared to the baselines.

Specifically, in sparse regression, our method using symmetry invariants is only slightly better than PySINDy in the reaction-diffusion system, but constantly succeeds in the Boussinesq equation where PySINDy fails. The failure of PySINDy is because the $u_{x}^{2}$ term in ([7](https://arxiv.org/html/2505.12083v1#S4.E7)) is not supported by its function library, showing that SINDy’s success relies heavily on the choice of function library. On the other hand, by enforcing the equation to be expressed in invariants, our method automatically identifies the proper function library.
[Section D.1](https://arxiv.org/html/2505.12083v1#A4.SS1) provides results for other variants of sparse regression.

For GP-based methods, [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1) displays the results with a fixed number of GP iterations for each dataset. We also include results with different numbers of iterations in [Section D.2](https://arxiv.org/html/2505.12083v1#A4.SS2). Generally, GP with invariants can identify the correct equation with fewer iterations and is considered more efficient.

### 4.4 Results on Noisy Data and Imperfect Symmetry

Figure: Figure 3: 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 SINDy under the same number of test functions. For systems with imperfect symmetry, strictly enforcing symmetry (SI) can hurt performance, but a relaxed symmetry constraint (SI-relaxed) is still better than no inductive bias (SINDy).

We test the robustness of our method under two challenging scenarios: (1) noise in observed data, and (2) PDE with imperfect symmetry.

In the first experiment, we add different levels of white noise to the simulated solution of the reaction-diffusion system. Since the derivatives estimated by finite difference is inaccurate with the noisy solution, we use the weak formulation of SINDy (Messenger & Bortz, 2021a), which does not require derivative estimation. The success probabilities of our method (SI) and SINDy are shown in [Figure 3](https://arxiv.org/html/2505.12083v1#S4.F3) (left), where $K$ is the number of test functions in weak SINDy. With the same $K$, our method consistently achieves higher success probability at different noise levels. Notably, when the noise level is high, our symmetry-constrained model performs better with fewer test functions ($K=100$).

In the second experiment, we simulate the two variants of ([10](https://arxiv.org/html/2505.12083v1#S4.E10)) (unequal diffusivities and external forcing) with different values for the symmetry-breaking parameter $\epsilon$ and add $2\%$ noise to the numerical solutions.
We compare three models: (1) our model with strictly enforced symmetry (SI), (2) our model with relaxed symmetry (SI-relaxed), and (3) weak SINDy as the baseline. The results for the two systems with symmetry breaking are shown in [Figure 3](https://arxiv.org/html/2505.12083v1#S4.F3) (center & right). As expected, SI has a much lower success probability when the symmetry-breaking factor becomes significant. Meanwhile, SI-relaxed remains highly competitive. It also has a clear 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 a model without any knowledge of symmetry.

More comprehensive results, e.g., samples of discovered equations, are provided in [Appendix D](https://arxiv.org/html/2505.12083v1#A4).

## 5 Discussion

In this paper, we propose to enforce symmetry in general methods for discovering symbolic differential equations by using the differential invariants of the symmetry group as the variable set in symbolic regression algorithms. We implement this general strategy in different classes of algorithms and observe improved accuracy, efficiency and robustness 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, our current framework cannot be applied if symmetry is 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 et al., 2024), which we leave for future work.

Another caveat of our method is the calculation of differential invariants. While solving for $\mathbf{v}^{(n)}(\eta)=0$ and applying the formula ([6](https://arxiv.org/html/2505.12083v1#S3.E6)) 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 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4)). Fortunately, this only requires a one-time effort. Once we have derived the invariants for a symmetry group, the results can be reused for any equation admitting the same symmetry.

## References

- Akhound-Sadegh et al. (2023)
Akhound-Sadegh, T., Perreault-Levasseur, L., Brandstetter, J., Welling, M., and Ravanbakhsh, S.
Lie point symmetry and physics informed networks.
*arXiv preprint arXiv:2311.04293*, 2023.
- Bakarji et al. (2022)
Bakarji, J., Callaham, J., Brunton, S. L., and Kutz, J. N.
Dimensionally consistent learning with buckingham pi.
*Nature Computational Science*, 2:834–844, 12 2022.
ISSN 2662-8457.
doi: 10.1038/s43588-022-00355-5.
- Biggio et al. (2021)
Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., and Parascandolo, G.
Neural symbolic regression that scales.
In Meila, M. and Zhang, T. (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.
- Brandstetter et al. (2022)
Brandstetter, J., Welling, M., and Worrall, D. E.
Lie point symmetry data augmentation for neural pde solvers.
In *International Conference on Machine Learning*, pp.  2241–2256. PMLR, 2022.
- Brunton et al. (2016)
Brunton, S. L., Proctor, J. L., and Kutz, J. N.
Discovering governing equations from data by sparse identification of nonlinear dynamical systems.
*Proceedings of the National Academy of Sciences*, 113(15):3932–3937, 2016.
doi: 10.1073/pnas.1517384113.
- Champion et al. (2019)
Champion, K., Lusch, B., Kutz, J. N., and Brunton, S. L.
Data-driven discovery of coordinates and governing equations.
*Proceedings of the National Academy of Sciences*, 116(45):22445–22451, 2019.
doi: 10.1073/pnas.1906995116.
- Cranmer (2023)
Cranmer, M.
Interpretable machine learning for science with pysr and symbolicregression.jl, 2023.
- Cranmer et al. (2020)
Cranmer, M., Sanchez Gonzalez, A., Battaglia, P., Xu, R., Cranmer, K., Spergel, D., and Ho, S.
Discovering symbolic models from deep learning with inductive biases.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), *Advances in Neural Information Processing Systems*, volume 33, pp.  17429–17442. Curran Associates, Inc., 2020.
- Cranmer et al. (2019)
Cranmer, M. D., Xu, R., Battaglia, P., and Ho, S.
Learning symbolic physics with graph networks, 2019.
- Dalton et al. (2024)
Dalton, D., Husmeier, D., and Gao, H.
Physics and lie symmetry informed gaussian processes.
In *Forty-first International Conference on Machine Learning*, 2024.
- de Silva et al. (2020)
de Silva, B., Champion, K., Quade, M., Loiseau, J.-C., Kutz, J., and Brunton, S.
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.21105/joss.02104.
URL [https://doi.org/10.21105/joss.02104](https://doi.org/10.21105/joss.02104).
- Dubčáková (2011)
Dubčáková, R.
Eureqa: software review, 2011.
- Finzi et al. (2021)
Finzi, M., Benton, G., and Wilson, A. G.
Residual pathway priors for soft equivariance constraints.
In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), *Advances in Neural Information Processing Systems*, volume 34, pp.  30037–30049. Curran Associates, Inc., 2021.
URL [https://proceedings.neurips.cc/paper_files/paper/2021/file/fc394e9935fbd62c8aedc372464e1965-Paper.pdf](https://proceedings.neurips.cc/paper_files/paper/2021/file/fc394e9935fbd62c8aedc372464e1965-Paper.pdf).
- Gaucel et al. (2014)
Gaucel, S., Keijzer, M., Lutton, E., and Tonda, A.
Learning dynamical systems using standard symbolic regression.
In Nicolau, M., Krawiec, K., Heywood, M. I., Castelli, M., García-Sánchez, P., Merelo, J. J., Rivas Santos, V. M., and Sim, K. (eds.), *Genetic Programming*, pp.  25–36, Berlin, Heidelberg, 2014. Springer Berlin Heidelberg.
- Grayeli et al. (2024)
Grayeli, A., Sehgal, A., Costilla Reyes, O., Cranmer, M., and Chaudhuri, S.
Symbolic regression with a learned concept library.
*Advances in Neural Information Processing Systems*, 37:44678–44709, 2024.
- Grundner et al. (2023)
Grundner, A., Beucler, T., Gentine, P., and Eyring, V.
Data-driven equation discovery of a cloud cover parameterization.
*arXiv preprint arXiv:2304.08063*, 2023.
- Holt et al. (2023)
Holt, S., Qian, Z., and van der Schaar, M.
Deep generative symbolic regression.
*arXiv preprint arXiv:2401.00282*, 2023.
- Kaheman et al. (2020)
Kaheman, K., Kutz, J. N., and Brunton, S. L.
Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics.
*Proceedings of the Royal Society A*, 476(2242):20200279, 2020.
- Kamienny et al. (2022)
Kamienny, P.-A., d’Ascoli, S., Lample, G., and Charton, F.
End-to-end symbolic regression with transformers.
*Advances in Neural Information Processing Systems*, 35:10269–10281, 2022.
- Kaptanoglu et al. (2022)
Kaptanoglu, A. A., de Silva, B. M., Fasel, U., Kaheman, K., Goldschmidt, A. J., Callaham, J., Delahunt, C. B., Nicolaou, Z. G., Champion, K., Loiseau, J.-C., Kutz, J. N., and Brunton, S. L.
Pysindy: A comprehensive python package for robust sparse system identification.
*Journal of Open Source Software*, 7(69):3994, 2022.
doi: 10.21105/joss.03994.
URL [https://doi.org/10.21105/joss.03994](https://doi.org/10.21105/joss.03994).
- Ko et al. (2024)
Ko, G., Kim, H., and Lee, J.
Learning infinitesimal generators of continuous symmetries from data.
*arXiv preprint arXiv:2410.21853*, 2024.
- Lee et al. (2022)
Lee, K., Trask, N., and Stinis, P.
Structure-preserving sparse identification of nonlinear dynamics for data-driven modeling.
In Dong, B., Li, Q., Wang, L., and Xu, Z.-Q. J. (eds.), *Proceedings of Mathematical and Scientific Machine Learning*, volume 190 of *Proceedings of Machine Learning Research*, pp.  65–80. PMLR, 15–17 Aug 2022.
- Martius & Lampert (2016)
Martius, G. and Lampert, C. H.
Extrapolation and learning equations.
*arXiv preprint arXiv:1610.02995*, 2016.
- Merler et al. (2024)
Merler, M., Haitsiukevich, K., Dainese, N., and Marttinen, P.
In-context symbolic regression: Leveraging large language models for function discovery.
*arXiv preprint arXiv:2404.19094*, 2024.
- Messenger & Bortz (2021a)
Messenger, D. A. and Bortz, D. M.
Weak sindy for partial differential equations.
*Journal of Computational Physics*, 443:110525, 2021a.
- Messenger & Bortz (2021b)
Messenger, D. A. and Bortz, D. M.
Weak sindy: Galerkin-based data-driven model selection.
*Multiscale Modeling & Simulation*, 19(3):1474–1497, 2021b.
- Mialon et al. (2023)
Mialon, G., Garrido, Q., Lawrence, H., Rehman, D., LeCun, Y., and Kiani, B.
Self-supervised learning with lie symmetries for partial differential equations.
*Advances in Neural Information Processing Systems*, 36:28973–29004, 2023.
- Newell (1985)
Newell, A. C.
*Solitons in mathematics and physics*.
SIAM, 1985.
- Olver (1993)
Olver, P. J.
*Applications of Lie groups to differential equations*, volume 107.
Springer Science & Business Media, 1993.
- Olver (1995)
Olver, P. J.
*Equivalence, invariants and symmetry*.
Cambridge University Press, 1995.
- Otto et al. (2023)
Otto, S. E., Zolman, N., Kutz, J. N., and Brunton, S. L.
A unified framework to enforce, discover, and promote symmetry in machine learning, 2023.
- Petersen et al. (2019)
Petersen, B. K., Landajuela, M., Mundhenk, T. N., Santiago, C. P., Kim, S. K., and Kim, J. T.
Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.
*arXiv preprint arXiv:1912.04871*, 2019.
- Qian et al. (2022)
Qian, Z., Kacprzyk, K., and van der Schaar, M.
D-code: Discovering closed-form odes from observed trajectories.
In *International Conference on Learning Representations*, 2022.
- Rao et al. (2022)
Rao, C., Ren, P., Liu, Y., and Sun, H.
Discovering nonlinear pdes from scarce data with physics-encoded learning.
*arXiv preprint arXiv:2201.12354*, 2022.
- Rudy et al. (2017)
Rudy, S. H., Brunton, S. L., Proctor, J. L., and Kutz, J. N.
Data-driven discovery of partial differential equations.
*Science advances*, 3(4):e1602614, 2017.
- Sahoo et al. (2018)
Sahoo, S., Lampert, C., and Martius, G.
Learning equations for extrapolation and control.
In *International Conference on Machine Learning*, pp.  4442–4450. PMLR, 2018.
- Schmidt & Lipson (2009)
Schmidt, M. and Lipson, H.
Distilling free-form natural laws from experimental data.
*science*, 324(5923):81–85, 2009.
- Shojaee et al. (2024)
Shojaee, P., Meidani, K., Gupta, S., Farimani, A. B., and Reddy, C. K.
Llm-sr: Scientific equation discovery via programming with large language models.
*arXiv preprint arXiv:2404.18400*, 2024.
- Shojaee et al. (2025)
Shojaee, P., Nguyen, N.-H., Meidani, K., Farimani, A. B., Doan, K. D., and Reddy, C. K.
Llm-srbench: A new benchmark for scientific equation discovery with large language models.
*arXiv preprint arXiv:2504.10415*, 2025.
- Takamoto et al. (2022)
Takamoto, M., Praditia, T., Leiteritz, R., MacKinlay, D., Alesiani, F., Pflüger, D., and Niepert, M.
Pdebench: An extensive benchmark for scientific machine learning.
*Advances in Neural Information Processing Systems*, 35:1596–1611, 2022.
- Udrescu & Tegmark (2020)
Udrescu, S.-M. and Tegmark, M.
Ai feynman: A physics-inspired method for symbolic regression.
*Science Advances*, 6(16):eaay2631, 2020.
- Udrescu et al. (2020)
Udrescu, S.-M., Tan, A., Feng, J., Neto, O., Wu, T., and Tegmark, M.
Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity.
*Advances in Neural Information Processing Systems*, 33:4860–4871, 2020.
- Wang et al. (2021)
Wang, R., Walters, R., and Yu, R.
Incorporating symmetry into deep dynamics models for improved generalization.
In *International Conference on Learning Representations*, 2021.
- Wang et al. (2022)
Wang, R., Walters, R., and Yu, R.
Approximately equivariant networks for imperfectly symmetric dynamics.
In *International Conference on Machine Learning*. PMLR, 2022.
- Wang et al. (2019)
Wang, Y., Wagner, N., and Rondinelli, J. M.
Symbolic regression in materials science.
*MRS Communications*, 9(3):793–805, 2019.
- Xie et al. (2022)
Xie, X., Samaei, A., Guo, J., Liu, W. K., and Gan, Z.
Data-driven discovery of dimensionless numbers and governing laws from scarce measurements.
*Nature Communications*, 13(1):7562, 2022.
doi: 10.1038/s41467-022-35084-w.
- Yang et al. (2024)
Yang, J., Rao, W., Dehmamy, N., Walters, R., and Yu, R.
Symmetry-informed governing equation discovery.
In *Advances in Neural Information Processing Systems (NeurIPS)*, 2024.
- Zhang et al. (2023)
Zhang, Z.-Y., Zhang, H., Zhang, L.-S., and Guo, L.-L.
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.

## Appendix A Related Works

#### Symbolic Regression.

Given the dataset $\{(x^{i},y^{i})\}\subset X\times Y$, symbolic regression (SR) aims to model the function $y=f(x)$ by a 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 modeling (Grundner et al., 2023), cosmology (Cranmer et al., 2020), etc. Various software have been developed for GP-based symbolic regression, e.g. Eureqa (Dubčáková, 2011) and PySR (Cranmer, 2023).

Another class of methods is sparse regression (Brunton et al., 2016), which assumes the function to be discovered can be written as a linear combination of predefined candidate functions and solves for the coefficient matrix. It has also been extended to discover more general equations, such as 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 (2016); Sahoo et al. (2018) represents a few earliest 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 et 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. (2023) pre-train an encoder-decoder 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, AI Feynman (Udrescu & Tegmark, 2020; Udrescu et al., 2020) uses properties like separability and compositionality to simplify the data. Cranmer et al. (2020) specifies 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 search. 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, evidenced by a recent benchmark specially designed for LLM-SR (Shojaee et 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 directly 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. A 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 et al., 2022; Lee et al., 2022).

#### 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 (Akhound-Sadegh et al., 2023; Zhang et al., 2023; Dalton et al., 2024), and self-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). 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). 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.

## Appendix B Math

### B.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 miscellaneous.**
| Symbols | Descriptions |
| --- | --- |
| $p$ | Number of independent variables of a PDE. |
| $q$ | Number of dependent variables of a PDE. |
| $X$ | Space of independent variables of a PDE: $X\subset\mathbb{R}^{p}$. Also used to denote the feature space of SR algorithms. |
| $U$ | Space of dependent variables of a PDE: $U\subset\mathbb{R}^{q}$. Assumed to be 1-dimensional unless otherwise stated. |
| $E$ | Total space of all variables of a PDE: $E=X\times U$. |
| $U_{k}$ | Space of strictly $k$th-order partial derivatives of variables in $U$ w.r.t variables in $X$. |
| $U^{(n)}$ | Space of all partial derivatives up to $n$th order (including the original variables in $U$): $U^{(n)}=U\times U_{1}\times\cdots\times U_{n}$. |
| $M^{(n)}$ | $n$th-order jet space: $M^{(n)}\subset X\times U^{(n)}$. |
| $T\mathcal{M}$ | The tangent bundle of a manifold $\mathcal{M}$. |
| $\mathbf{x}$ | Independent variables of a PDE: $\mathbf{x}\in\mathbb{R}^{p}$. |
| $t$ | Time variable. |
| $x,y$ | Spatial variables in PDE contexts. Also used to denote the features and labels of SR algorithms, where $x$ can denote multi-dimensional features. |
| $u,\mathbf{u}$ | Dependent variable(s) of a PDE: $u\in\mathbb{R}$ and $\mathbf{u}\in\mathbb{R}^{q}$. |
| $u^{(n)},\mathbf{u}^{(n)}$ | The collection of all up to $n$-th order partial derivatives of $u$ or $\mathbf{u}$. |
| $df$ | The (ordinary) differential of a function. For a differential function $f:M^{(n)}\to\mathbb{R}$, $df=\sum_{j}\frac{\partial f}{\partial x^{j}}dx^{j}+\sum_{\alpha}\frac{\partial<br>f<br>}{\partial u_{\alpha}}du_{\alpha}$. |
| $D_{i}f$ | The total derivative of a differential function $f:M^{(n)}\to\mathbb{R}$ w.r.t the $i$th independent variable. For example, if $p=q=1$, $D_{1}f=\frac{\partial f}{\partial x}+\sum_{k=0}^{\infty}u_{k+1}\frac{\partial f<br>}{\partial u_{k}}$, where $u_{k}\coloneqq\partial^{k}u/\partial x^{k}$. |
| $Df$ | The total differential of a differential function $f:M^{(n)}\to\mathbb{R}$, i.e. $Df=D_{i}f\ dx^{i}$. |
| $g$ | A group element with an action on $E$ ([2](https://arxiv.org/html/2505.12083v1#S2.E2)). |
| $\mathbf{v}$ | A vector field on the total space $E$ ([3](https://arxiv.org/html/2505.12083v1#S2.E3)), representing an infinitesimal transformation. A list of multiple vector fields are indexed by subscripts. |
| $\mathrm{pr}^{(n)}g$ | $n$th-order prolongation of $g$ acting on $M^{(n)}$. |
| $\mathrm{pr}^{(n)}\mathbf{v}$ | $n$th-order prolongation of $\mathbf{v}$ acting on $M^{(n)}$. |
| $g^{(n)},\mathbf{v}^{(n)}$ | Equivalent to $\mathrm{pr}^{(n)}g$ and $\mathrm{pr}^{(n)}\mathbf{v}$, respectively. |
| $\mathrm{pr}\ \mathbf{v}$ | The (infinite) prolongation of $\mathbf{v}$. For an $n$th-order differential function $f(\mathbf{x},\mathbf{u}^{(n)})$, $\mathrm{pr}\ \mathbf{v}(f)=\mathrm{pr}^{(n)}\mathbf{v}(f)$. |
| $\eta,\zeta,\vartheta$ | Differential invariants of a symmetry group. $\eta$ is used by default. The other letters are used to distinguish between invariants of different orders. |
| $\ell,\bm{\ell}$ | The LHS of SINDy equation ([4](https://arxiv.org/html/2505.12083v1#S2.E4)). Often assumed to be time derivatives. |
| $\bm{\theta}$ | A column vector containing all SINDy library functions: $\bm{\theta}=[\theta^{1},\cdots,\theta^{m}]$ |
| $\mathbf{w},W$ | The SINDy parameters. For only one equation, $\mathbf{w}=[w^{1},\cdots,w^{m}]$ is a row vector. For multiple equations, $W=[w^{ij}]$ is a $q\times m$ matrix. |
| $\mathbf{X},\mathbf{y}$ | Concatenated matrix/vector of features/labels of all datapoints for symbolic regression. |
| $[N]$ | List of positive integers up to $N$, i.e. $[1,2,\cdots,N]$ for any $N\in\mathbb{Z}^{+}$. |
| $1:N$ | Equivalent to $[N]$. |
| LHS, RHS | Left- and Right-hand side of an equation. |

### B.2 Extended Background on PDE Symmetry

References for the below material include Olver (1993), Olver (1995).

#### Prolonged group actions

Let $E=X\times U\simeq\mathbb{R}^{p}\times\mathbb{R}^{q}$ be endowed with the action of a group $G$ via point transformations. Then group elements $g\in G$ act locally on functions $\mathbf{u}=f(\mathbf{x})$, therefore also on derivatives of these functions. This in turn induces, at least pointwise, “prolonged" transformations on jet spaces: $(\tilde{\mathbf{x}},\tilde{\mathbf{u}}^{(n)})=\mathrm{pr}^{(n)}g\cdot(\mathbf{
x},\mathbf{u}^{(n)})$.

Let $J=(j_{1},\dots,j_{n}),1\leq j_{\nu}\leq p$ be an $n$-tuple of indices of independent variables and $1\leq\alpha\leq q$. We will use the shorthand

$$ $u^{\alpha}_{J}:=\frac{\partial^{J}u^{\alpha}}{\partial x^{J}}:=\frac{\partial^ {|J|}u^{\alpha}}{\partial x_{j_{1}}\cdots\partial x_{j_{n}}}$ $$

and

$$ $D_{J}:=D_{j_{1}}\cdots D_{j_{n}}.$ $$

It is 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

$$ $\mathbf{v}=\sum_{i=1}^{p}\xi^{i}(\mathbf{x},\mathbf{u})\frac{\partial}{ \partial x^{i}}+\sum_{\alpha=1}^{q}\varphi^{\alpha}(\mathbf{x},\mathbf{u}) \frac{\partial}{\partial u^{\alpha}},$ $$

its characteristic is a $q$-tuple $Q=(Q^{1},\dots,Q^{q})$ of functions with

$$ $Q^{\alpha}(\mathbf{x},\mathbf{u}^{(1)})=\varphi^{\alpha}(\mathbf{x},\mathbf{u} )-\sum_{i=1}^{p}\xi^{i}(\mathbf{x},\mathbf{u})\frac{\partial u^{\alpha}}{ \partial x^{i}}.$ $$

Now the prolongation of $\mathbf{v}$ to order $n$ is defined by

$$ $\mathrm{pr}^{(n)}\mathbf{v}=\sum_{i=1}^{p}\xi^{i}(\mathbf{x},\mathbf{u})\frac{ \partial}{\partial x^{i}}+\sum_{\alpha=1}^{q}\sum_{\#J=n}\varphi_{J}^{\alpha}( \mathbf{x},\mathbf{u}^{(n)})\frac{\partial}{\partial u^{\alpha}_{J}}.$ (11) $$

Here $J$ ranges over all $n$-tuples $J=(j_{1},\dots,j_{n}),1\leq j_{\nu}\leq p$ and the $\varphi^{\alpha}_{J}$ are given by

$$ $\varphi^{\alpha}_{J}=D_{J}Q^{\alpha}+\sum_{i=1}^{p}\xi^{i}\mathbf{u}^{\alpha}_ {J,i}.$ $$

We remark that the prolongation of $\mathbf{v}$ has been described explicitly in terms of the coefficients of $\mathbf{v}$ and their derivatives.

### B.3 Proof of Proposition 3.3

Olver (1995) provides the following general theorem to construct higher-order differential invariants from a contact-invariant coframe. We refer the readers to Chapter 5 of Olver (1995) for definitions of relevant concepts, e.g., contact forms and contact-invariant forms and coframes.

###### Theorem B.1 (Thm. 5.48, (Olver, 1995 ) ) .

Let $G$ be a transformation group acting on a space with $p$ independent variables and $q$ dependent variables. Suppose $\omega^{1},...,\omega^{p}$ is a contact-invariant coframe for $G$, and let $\mathcal{D}_{j}$ be the associated invariant differential operators defined via $Df=D_{j}f\ dx^{j}=\mathcal{D}_{j}f\ \omega^{j}$. If there are a maximal number of independent, strictly $n$th-order differential invariants $\zeta^{1},\cdots,\zeta^{q_{n}}$, $q_{n}=\binom{p+n-1}{n}$, then the set of differentiated invariants $\mathcal{D}_{i}\zeta^{\nu}$, $i\in[p]$, $\nu\in[q_{n}]$, contains a complete set of independent, strictly $(n+1)$th-order differential invariants.

Specifically, the condition that there exist a maximal number of differential invariants of order exactly $n$ is guaranteed if $n$ is at least $\mathrm{dim}G$.

Our proposition is a derived result from the above theorem, which provides a concrete way of computation from lower-order invariants to higher-order ones:

###### Proposition B.2 .

Let $G$ be a local group acting on $X\times U\simeq\mathbb{R}^{p}\times\mathbb{R}$. Let $\eta^{1},\eta^{2},\cdots,\eta^{p}$ be any $p$ differential invariants of $G$ whose horizontal Jacobian $J=[D_{i}\eta^{j}]$ is non-degenerate on an open subset $\Omega\subset M^{(n)}$. If there are a maximal number of independent, strictly $n$th-order differential invariants $\zeta^{1},\cdots,\zeta^{q_{n}}$, $q_{n}=\binom{p+n-1}{n}$, then the following set contains a complete set of independent, strictly $(n+1)$th-order differential invariants defined on $\Omega$:

$$ $\frac{\mathrm{det}(D_{i}\tilde{\eta}^{j}_{(k,k^{\prime})})}{\mathrm{det}(D_{i} \eta^{j})},\ \forall k\in[p],k^{\prime}\in[q_{n}],$ (12) $$

where $i,j\in[p]$ are matrix indices, $D_{i}$ denotes the total derivative w.r.t $i$-th independent variable and $\tilde{\eta}^{j}_{(k,k^{\prime})}=[\eta^{1},...,\eta^{k-1},\zeta^{k^{\prime}},
\eta^{k+1},...,\eta^{p}]$.

###### Proof.

We show that the total differentials of the differential invariants $\eta^{1},...,\eta^{p}$ can be used to construct a contact-invariant coframe of $G$ and then derive the associated invariant differential operators to complete the proof.

First, note that for any differential invariant $\eta$ of $G$, its total differential $\omega=D\eta=D_{j}\eta\ dx^{j}$ can be written as

$$ $\omega=\omega_{o}+\theta,$ (13) $$

where $\omega_{o}\coloneqq d\eta=\sum_{i\in[p]}\frac{\partial F}{\partial x^{i}}dx^{i
}+\sum_{|\alpha|\leq n}\frac{\partial F}{\partial u_{\alpha}}du_{\alpha}$ is the ordinary differential of $\eta:M^{(n)}\to\mathbb{R}$ and $\theta$ is a contact form.

Since $\eta$ is a differential invariant, its differential $\omega_{o}=d\eta$ is an invariant one-form on $M^{(n)}$, i.e. $(g^{(n)})^{*}\omega_{o}=\omega_{o}$.

Also, a prolonged group action maps contact forms to contact forms. To see this, note that a prolonged group action $g^{(n)}$ maps the prolonged graph of any function to the prolonged graph of a transformed function. Then, for any contact form $\theta$, $(g^{(n)})^{*}\theta$ is annihilated by all prolonged functions $f^{(n)}$, thus a contact form by definition:

$$ $\displaystyle(f^{(n)})^{*}((g^{(n)})^{*}\theta)$ $\displaystyle=(g^{(n)}\circ f^{(n)})^{*}\theta$ (14) $\displaystyle=((g\cdot f)^{(n)})^{*}\theta$ (15) $\displaystyle=0.$ (16) $$

Then, from ([13](https://arxiv.org/html/2505.12083v1#A2.E13)), we have

$$ $\displaystyle(g^{(n+1)})^{*}\omega$ $\displaystyle=(g^{(n)})^{*}\omega_{o}+(g^{(n+1)})^{*}\theta$ (17) $\displaystyle=\omega_{o}+\theta^{\prime}$ (18) $\displaystyle=\omega+(\theta^{\prime}-\theta)$ (19) $$

where $\theta^{\prime}$ is some contact form and so is $\theta^{\prime}-\theta$. Thus, $\omega$ is contact-invariant. For the $p$ differential invariants $\eta^{1},\cdots,\eta^{p}$, we have $p$ contact-invariant one-forms $\omega^{1},\cdots,\omega^{p}$, respectively.

Next, we prove that $\omega^{1},\cdots,\omega^{p}$ are linearly independent and form a coframe. Assume there exists smooth coefficients $c^{j}$ such that $\sum_{j}c^{j}\omega^{j}=0$. Then, regrouping the coefficients of the horizontal forms $dx^{i}$, we have

$$ $0=\sum_{i,j}c^{j}D_{i}\eta^{j}dx^{i}=\sum_{i}\left(\sum_{j}c^{j}D_{i}\eta^{j} \right)dx^{i}.$ (20) $$

Because the $dx^{i}$ are linearly independent, each coefficient of $dx^{i}$ must vanish, i.e. $\tensor{J}{{}_{i}^{j}}c^{j}=0$. Since the Jacobian $J=[D_{i}\eta^{j}]$ is non-degenerate, the only solution is $c^{j}=0$ (on the open subset $\Omega\in M^{(n)}$). Thus, $\omega^{1},\cdots,\omega^{p}$ form a contact-invariant coframe. According to [Theorem B.1](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem1), the associated invariant differential operators of the coframe take a complete set of same-order invariants to a complete set of one-order-higher invariants.

The remaining step is to obtain the invariant differential operators explicitly in terms of $\eta^{j}$. Recall the formula in [Theorem B.1](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem1) that defines the invariant differential operators:

$$ $D_{i}f\ dx^{i}=\mathcal{D}_{j}f\ \omega^{j}.$ (21) $$

Expanding $\omega^{j}=D\eta^{j}=D_{i}\eta^{j}\ dx^{i}$, we have the following linear system of invariant differential operators $\mathcal{D}_{j}$:

$$ $\begin{bmatrix}D_{1}\\ D_{2}\\ \vdots\\ D_{p}\end{bmatrix}=\begin{bmatrix}D_{1}\eta^{1}&D_{1}\eta^{2}&\cdots&D_{1}\eta ^{p}\\ D_{2}\eta^{1}&D_{2}\eta^{2}&\cdots&D_{2}\eta^{p}\\ \vdots&\vdots&&\vdots\\ D_{p}\eta^{1}&D_{p}\eta^{2}&\cdots&D_{p}\eta^{p}\end{bmatrix}\begin{bmatrix} \mathcal{D}_{1}\\ \mathcal{D}_{2}\\ \vdots\\ \mathcal{D}_{p}\end{bmatrix}.$ (22) $$

Since $J=[D_{i}\eta^{j}]$ is non-degenerate, Cramer’s rule yields

$$ $\mathcal{D}_{k}\zeta=\frac{\mathrm{det}(D_{i}\eta^{1}\ |\ \cdots\ |\ D_{i}\eta ^{k-1}\ |\ D_{i}\zeta\ |\ D_{i}\eta^{k+1}\ |\ \cdots\ |\ D_{i}\eta^{p})}{ \mathrm{det}(D_{i}\eta^{j})}.$ (23) $$

∎

###### Remark B.3 .

We require that the differential invariants $\eta^{1},\cdots,\eta^{p}$ has a nondegenerate horizontal Jacobian $[D_{i}\eta^{j}]$, which is a 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 $[\partial_{i}\eta^{j}]$, where $i\in[q_{n}]$ indexes the jet space variables $(\mathbf{x},u^{(n)})$, but has a lower-rank horizontal Jacobian. For example, consider $\eta^{1}=u_{x}$ and $\eta^{2}=u_{y}$. In the full Jacobian, $\partial\eta^{j}/\partial u_{x}$ and $\partial\eta^{j}/\partial u_{y}$ form the identity, so it has full rank. However, its horizontal Jacobian containing total derivatives is given by $\begin{bmatrix}u_{xx}&u_{xy}\\
u_{xy}&u_{yy}\end{bmatrix}$, which is not invertible on the subset of the jet space where $u_{xx}u_{yy}-u_{xy}^{2}=0$.

In practice, this non-degeneracy condition can be easily checked once we have the symbolic expressions of the $p$ differential invariants.

###### Remark B.4 .

When $p=1$, [Proposition B.2](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem2) is equivalent to the following (Prop. 2.53, Olver (1993)):

If $y=\eta(x,u^{(n)})$ and $w=\zeta(x,u^{(n)})$ are $n$-th order differential invariants of $G$, then $\frac{dw}{dy}\equiv\frac{D_{x}\zeta}{D_{x}\eta}$ is an $(n+1)$-th order differential invariant of $G$.
Specifically, if $y=\eta(x,u)$ and $w=\zeta(x,u,u_{x})$ form a complete set of functionally independent differential invariants of $\mathrm{pr}^{(1)}G$, the complete set of functionally independent differential invariants for $\mathrm{pr}^{(n)}G$ is then given by

$$ $y,w,dw/dy,...,d^{n-1}w/dy^{n-1}.$ (24) $$

###### Theorem B.1 (Thm. 5.48, (Olver, 1995 ) ) .

###### Proposition B.2 .

###### Proof.

###### Remark B.3 .

###### Remark B.4 .

### B.4 Examples of Computing Differential Invariants

###### Example B.5 .

Consider the group $\mathrm{SO}(2)$ acting on $X\times U\simeq\mathbb{R}^{2}\times\mathbb{R}$ by standard rotation in the 2D space of independent variable and trivial action on $U$, i.e. its infinitesimal generator given by $\mathbf{v}=y\partial_{x}-x\partial_{y}$.

First, we solve for a complete set of the ordinary and first-order invariants. The two ordinary invariants are given by $\eta_{1}(x,y,u)=\frac{1}{2}(x^{2}+y^{2})$ and $\eta_{2}(x,y,u)=u$. ([6](https://arxiv.org/html/2505.12083v1#S3.E6)) dictates how we construct higher-order invariants using these two functionally independent invariants and another arbitrary invariant. For notational convenience, we convert ([6](https://arxiv.org/html/2505.12083v1#S3.E6)) to operators defined according to $\eta_{2}$ and $\eta_{1}$, respectively:

$$ $\displaystyle\mathcal{O}_{1}$ $\displaystyle=\frac{xD_{y}-yD_{x}}{xu_{y}-yu_{x}}$ (25) $\displaystyle\mathcal{O}_{2}$ $\displaystyle=\frac{u_{y}D_{x}-u_{x}D_{y}}{xu_{y}-yu_{x}}$ (26) $$

Then, we need to find another new differential invariant, because applying these operators on $\eta_{1}$ and $\eta_{2}$ leads to trivial results. Since $\eta_{1}$ and $\eta_{2}$ 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

$$ $\mathrm{pr}^{(1)}\mathbf{v}=\mathbf{v}+u_{y}\partial_{u_{x}}-u_{x}\partial_{u_ {y}}$ (27) $$

Solving for $\mathrm{pr}^{(1)}\mathbf{v}$ gives two first-order invariants, $\zeta_{1}=xu_{y}-yu_{x}$ and $\zeta_{2}=xu_{x}+yu_{y}$. Note that the differential invariant $\zeta_{1}$ is exactly the common denominator in $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, so we can simplify $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ by using only their numerators, i.e.

$$ $\displaystyle\mathcal{O}_{1}$ $\displaystyle=xD_{y}-yD_{x}$ (28) $\displaystyle\mathcal{O}_{2}$ $\displaystyle=u_{y}D_{x}-u_{x}D_{y}$ (29) $$

Note that $\mathcal{O}_{2}$ has first-order coefficients, which may complicate things in the subsequent calculation. Denoting the space of all continuous functions of the existing four invariants as $\mathcal{I}=\mathcal{C}(\eta_{1},\eta_{2},\zeta_{1},\zeta_{2})$, we can choose any new operator within the $\mathcal{I}$-module spanned by $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ that makes things easier. Specifically, we use the following operator

$$ $\displaystyle\tilde{\mathcal{O}}_{2}$ $\displaystyle=\frac{\zeta_{2}}{\zeta_{1}}\mathcal{O}_{1}+\frac{2\eta_{1}}{ \zeta_{1}}\mathcal{O}_{2}$ (30) $\displaystyle=xD_{x}+yD_{y}$ (31) $$

Then, we apply these operators to the first-order invariants, which raise the order by one and give us the second-order invariants. For example, applying $O_{1}$ to $\zeta_{1}$, we have

$$ $\displaystyle\mathcal{O}_{1}\zeta_{1}$ $\displaystyle=xD_{y}\zeta_{1}-yD_{x}\zeta_{1}$ (32) $\displaystyle=x(xu_{yy}-u_{x}-yu_{xy})-y(u_{y}+xu_{xy}-yu_{xx})$ (33) $\displaystyle=x^{2}u_{yy}+y^{2}u_{xx}-xu_{x}-yu_{y}-2xyu_{xy}$ (34) $$

Note that $\zeta_{2}=xu_{x}+yu_{y}$ is a first-order invariant, so we can further remove it from the formula and get a simplified second-order invariant

$$ $\vartheta_{1}=x^{2}u_{yy}+y^{2}u_{xx}-2xyu_{xy}$ (35) $$

Similarly, we compute $\mathcal{O}_{1}\zeta_{2}$, $\tilde{\mathcal{O}}_{2}\zeta_{1}$ and $\tilde{\mathcal{O}}_{2}\zeta_{2}$ and obtain the following, respectively:

$$ $\displaystyle\vartheta_{2}=\vartheta_{3}$ $\displaystyle=\zeta_{1}+xy(u_{yy}-u_{xx})+(x^{2}-y^{2})u_{xy}$ (36) $\displaystyle\equiv xy(u_{yy}-u_{xx})+(x^{2}-y^{2})u_{xy}$ (37) $\displaystyle\vartheta_{4}$ $\displaystyle=\zeta_{2}+x^{2}u_{xx}+y^{2}u_{yy}+2xyu_{xy}$ (38) $\displaystyle\equiv x^{2}u_{xx}+y^{2}u_{yy}+2xyu_{xy}$ (39) $$

The above 8 invariants should form a complete set of second-order differential invariants of $\mathbf{v}=x\partial_{y}-y\partial_{x}$. To verify, note that the Laplacian $\Delta u=u_{xx}+u_{yy}$, which is a well-known rotational invariant, can be written in terms of these differential functions:

$$ $\displaystyle\Delta u=u_{xx}+u_{yy}$ $\displaystyle=\frac{(x^{2}+y^{2})(u_{xx}+u_{yy})}{x^{2}+y^{2}}$ (40) $\displaystyle=\frac{\vartheta_{1}+\vartheta_{4}}{2\eta_{1}}$ (41) $$

Another second-order rotational invariant, the trace of the squared Hessian matrix, $u_{xx}^{2}+2u_{xy}^{2}+u_{yy}^{2}$, is recovered by

$$ $u_{xx}^{2}+2u_{xy}^{2}+u_{yy}^{2}=\frac{\vartheta_{1}^{2}+2\vartheta_{2}^{2}+ \vartheta_{4}^{2}}{4\eta_{1}^{2}}$ (42) $$

On the other hand, these 8 invariants are apparently not functionally independent - note that $\vartheta_{2}=\mathcal{O}_{1}\zeta_{2}$ and $\vartheta_{3}=\tilde{\mathcal{O}}_{2}\zeta_{1}$ are the same. While this may be some coincidence, eventually it is not surprising because we would expect to see 3 functionally independent strictly second-order differential invariants instead of 4, since $(u_{xx},u_{yy},u_{xy})\in U_{2}$ is only 3-dimensional.

###### Example B.6 (Scaling and translation) .

Consider the vector field $\mathbf{v}_{1}=t\partial_{t}+ax\partial_{x}+bu\partial_{u}$. It generates the scaling symmetry $t\mapsto\lambda t,x\mapsto\lambda^{a}x,u\mapsto\lambda^{b}u$. The ordinary invariants of this symmetry are $t^{b}u^{-1}$ and $x^{a}u^{-1}$.
The higher-order invariants are given by
$\eta_{(\alpha,\beta)}=x^{\alpha}t^{\beta}u_{x^{(\alpha)}t^{(\beta)}}u^{-1}$,
where $\alpha$ and $\beta$ denote the orders of partial derivatives w.r.t $t$ and $x$, e.g. $u_{x^{(2)}t^{(1)}}\coloneqq u_{xxt}$.

Besides the scaling symmetry, we can consider other common symmetries simultaneously, e.g. translation symmetries in both space and time, $\mathbf{v}_{2}=\partial_{x}$ and $\mathbf{v}_{3}=\partial_{t}$. These symmetries, along with the scaling symmetry $\mathbf{v}_{1}$, span a 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

$$ $\eta_{(\alpha,\beta)}=u_{x^{(\alpha)}t^{(\beta)}}u_{x}^{\frac{b-a\alpha-\beta} {a-b}},\ \alpha\geq 0,\beta\geq 0.$ (43) $$

###### Example B.5 .

###### Example B.6 (Scaling and translation) .

### B.5 Proof of Proposition 3.4

[Proposition 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4), restated below, aligns our symmetry constraint into the SINDy framework and results in a set of constraints on the SINDy parameters.

###### Proposition B.7 .

Let $\bm{\ell}(\mathbf{x},\mathbf{u}^{(n)})=W\bm{\theta}(\mathbf{x},\mathbf{u}^{(n)})$ be a system of $q$ differential equations admitting a symmetry group $G$, where $\mathbf{x}\in\mathbb{R}^{p}$, $\mathbf{u}\in\mathbb{R}^{q}$, $\bm{\theta}\in\mathbb{R}^{m}$. Assume that there exist some $n$th-order invariants of $G$, $\eta_{0}^{1:q}\text{ and }\eta^{1:K}$, s.t. (1) the system of differential equations can be expressed as $\bm{\eta}_{0}=W^{\prime}\bm{\theta}^{\prime}(\bm{\eta})$, where $\bm{\eta}_{0}=[\eta_{0}^{1:q}]$ and $\bm{\eta}=[\eta^{1:K}]$, and (2) $\eta_{0}^{i}=T^{ijk}\theta^{k}\ell^{j}$ and $(\theta^{\prime})^{i}=S^{ij}\theta^{j}$, for some library functions $\bm{\theta}^{\prime}(\bm{\eta})$ and some constant tensors $W^{\prime}$, $T$ and $S$. Then, the space of all possible $W$ is a linear subspace of $\mathbb{R}^{q\times m}$.

###### 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

$$ $\ell^{i}=W^{ij}\theta^{j}.$ (44) $$

Combining the conditions about the differential invariants, we know that the equation can be equivalently expressed as

$$ $T^{ijk}\theta^{k}\ell^{j}=(W^{\prime})^{ij}S^{jk}\theta^{k}$ (45) $$

for some $W^{\prime}\in\mathbb{R}^{q\times m^{\prime}}$, where $m^{\prime}$ is the number of invariant functions in $\bm{\theta}^{\prime}$.

Substituting ([44](https://arxiv.org/html/2505.12083v1#A2.E44)) into ([45](https://arxiv.org/html/2505.12083v1#A2.E45)) and rearranging the indices, the principle of symmetry invariants then translates to the following constraint on $W$: there exists some $W^{\prime}\in\mathbb{R}^{q\times m^{\prime}}$ s.t.

$$ $\tensor{T}{{}_{i}^{r}{}^{k}}\theta_{k}\tensor{W}{{}_{r}^{l}}\theta_{l}=(W^{ \prime}\tensor{)}{{}_{i}^{k}}\tensor{S}{{}_{k}^{j}}\theta_{j},\forall\mathbf{x },\mathbf{u}^{(n)}.$ (46) $$

To solve for $W$, we first eliminate the dependency on the variables $\mathbf{x}$ and $\mathbf{u}^{(n)}$ from the equation. We adopt a procedure similar to Yang et al. (2024). Denote $\mathbf{z}=(\mathbf{x},\mathbf{u}^{(n)})$. Define a functional $M_{\bm{\theta}}$ as mapping a function to its coordinate in the function space spanned by ${\bm{\theta}}$, i.e. $M_{\bm{\theta}}:(\mathbf{z}\mapsto\tensor{c}{{}^{j}}\theta_{j}(\mathbf{z}))
\mapsto(c^{1},c^{2},\cdots,c^{m})$. Before we proceed, note that the LHS of ([46](https://arxiv.org/html/2505.12083v1#A2.E46)) contains the products of functions $\theta_{k}(\mathbf{z})\theta_{l}(\mathbf{z})$, which may or may not be included in the original function library $\bm{\theta}$. Therefore, we denote $\tilde{\bm{\theta}}(\mathbf{z})=[{\bm{\theta}}(\mathbf{z})\ ||\ \{\theta_{k}
\theta_{l}\notin{\bm{\theta}}\}]$ as the collection of all library functions $\theta_{k}$ and all their products $\theta_{k}\theta_{l}$. The invariant functions ${\bm{\theta}}^{\prime}(\bm{\eta})$ can also be rewritten in terms of the prolonged library: ${\bm{\theta}}^{\prime}(\bm{\eta})=\tilde{S}\tilde{\bm{\theta}}$, where $\tilde{S}_{1:m}=S$.

Then, applying $M_{\tilde{\bm{\theta}}}$ to ([46](https://arxiv.org/html/2505.12083v1#A2.E46)), we have

$$ $M_{\tilde{\bm{\theta}}}(\tensor{T}{{}_{i}^{r}{}^{k}}\theta_{k}\tensor{W}{{}_{r }^{l}}\theta_{l})=(W^{\prime}\tensor{)}{{}_{i}^{k}}\tensor{\tilde{S}}{{}_{k}^{ j}}.$ (47) $$

Further expanding the LHS, we have

$$ $\tensor{T}{{}_{i}^{r}{}^{k}}\tensor{W}{{}_{r}^{l}}\tensor{\Gamma}{{}_{k}{}_{l} ^{j}}=(W^{\prime}\tensor{)}{{}_{i}^{k}}\tensor{\tilde{S}}{{}_{k}^{j}},$ (48) $$

where $\Gamma$ satisfies $\theta_{k}\theta_{l}=\tensor{\Gamma}{{}_{k}{}_{l}^{j}}\tilde{\theta}_{j}$. In other words, the rows of the LHS fall in the row space of $\tilde{S}$. Let $\tilde{S}^{\perp}$ be the basis matrix for the null space of $\tilde{S}$, i.e. $\tilde{S}\tilde{S}^{\perp}=0$, we have

$$ $\tensor{T}{{}_{i}^{r}{}^{k}}\tensor{W}{{}_{r}^{l}}\tensor{\Gamma}{{}_{k}{}_{l} ^{j}}(\tilde{S}^{\perp}\tensor{)}{{}_{j}{}_{s}}=0,$ (49) $$

suggesting that $W$ must lie in a linear subspace of $\mathbb{R}^{q\times m}$.

∎

###### Remark B.8 .

In practice, to solve for ([49](https://arxiv.org/html/2505.12083v1#A2.E49)), we first rearrange ([49](https://arxiv.org/html/2505.12083v1#A2.E49)) into $M\mathrm{vec}(W)=0$, where $M$ has shape $(\tilde{S}\text{.shape}[2]\times q,\ q\times m)$. Then, we perform SVD on $M$ and apply a threshold of $10^{-6}$ to the singular values. The right singular vectors corresponding to the singular values smaller than the threshold then form a basis of the linear subspace $\mathrm{vec}(W)$ lies in.

###### Proposition B.7 .

###### Proof.

###### Remark B.8 .

## Appendix C Implementation Details

This section discusses some detailed considerations in implementing the sparse regression-based methods described in [Section 3.3](https://arxiv.org/html/2505.12083v1#S3.SS3) and [3.4](https://arxiv.org/html/2505.12083v1#S3.SS4). Contents include:

- •
[Section C.1](https://arxiv.org/html/2505.12083v1#A3.SS1): An algorithmic description of direct sparse regression with symmetry invariants.
- •
[Section C.2](https://arxiv.org/html/2505.12083v1#A3.SS2): Converting the symmetry invariant condition as linear constraints on the sparse regression parameters.
- •
[Section C.3](https://arxiv.org/html/2505.12083v1#A3.SS3): Using differential invariants in weak SINDy via the linear constraints, as well as other considerations.

### C.1 Direct Sparse Regression With Symmetry Invariants

The first approach to enforcing symmetry in sparse regression, as discussed in [Section 3.3](https://arxiv.org/html/2505.12083v1#S3.SS3), is to directly use the symmetry invariants as the variables and their functions specified by a function library as the RHS features. Similar to [Algorithm 1](https://arxiv.org/html/2505.12083v1#alg1) for general symbolic regression methods, we provide a detailed algorithm for sparse regression below. Following the setup from SINDy, we aim to discover a system of $q$ differential equations for $q$ dependent variables.

Figure: Algorithm 2 Sparse regression with symmetry invariants

The configuration from the original SINDy model, i.e., the LHS $\bm{\ell}$ and the function library $\{\theta^{j}\}$, are used to construct a new equation model in terms of the invariants.
It 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., 2022) implementation, a function $\theta$ is provided in a lambda format lambda x, y: x * y. Thus, $\theta$ can be applied to both the original variables, e.g. $\theta(z_{1},z_{2})=z_{1}z_{2}$, and the invariant functions, e.g. $\theta(\eta_{1},\eta_{2})=\eta_{1}\eta_{2}$.

### C.2 Symmetry Invariant Condition as Linear Constraints

Instead of directly using the invariant functions $\eta$ as the features and labels for regression, we can derive a set of linear constraints from the fact that the equation can be rewritten in terms of invariant functions. As shown in [Section B.5](https://arxiv.org/html/2505.12083v1#A2.SS5), a 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,q,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 $W^{jk}=Q^{ijk}\beta^{i}$, where $\beta$ is the learnable parameter, and discover the equation using the original SINDy objective as described in [Section 2.2](https://arxiv.org/html/2505.12083v1#S2.SS2).

In practice, we observe that the basis $Q$ obtained from SVD is not sparse. Indeed, SVD does not inherently encourage sparsity in the singular vectors. The lack of sparsity can pose a problem when we perform sequential thresholding in sparse regression. Specifically, in SINDy, the entries in $W$ that are close to zero are filtered out at the end of each iteration, which serves as a proxy to the $L_{0}$ regularization. Since we fix $Q$ and only optimize $\beta$, a straightforward modification to the sequential thresholding procedure is to threshold the entries in $\beta$ instead of those in $W$. However, if $Q$ is dense, even a sparse vector $\beta$ can lead to a dense $W$, which contradicts the purpose of sparse regression.

Therefore, after performing SVD, we apply a Sparse PCA to $Q$ to obtain a sparsified basis, also of shape $(r,q,m)$:

[Figure 4](https://arxiv.org/html/2505.12083v1#A3.F4) shows an example of the original basis solved from SVD (top $7\times 2$ grid) and the sparsified basis using sparse PCA (bottom $7\times 2$ grid). This is used in our experiment on the reaction-diffusion system ([10](https://arxiv.org/html/2505.12083v1#S4.E10)).

Figure: Figure 4: Basis for the SINDy parameter subspace that preserves $\mathrm{SO}(2)$ symmetry $\mathbf{v}=-v\partial_{u}+u\partial_{v}$. The SINDy parameter $W$ has dimension $2\times 19$. The two rows correspond to the two equations with $u_{t}$ and $v_{t}$ as the LHSs. The RHS contains $19$ features, including all monomials of $u,v$ 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,u^{2}+v^{2}\}\bigcup\{\mathbf{u}\cdot\mathbf{u}_{\mu}\}\bigcup\{ \mathbf{u}^{\perp}\cdot\mathbf{u}_{\mu}\}$, where $\mathbf{u}=(u,v)^{T}$ and $\mu$ is a multiindex of $t,x,y$ with order no more than $2$. The top $7\times 2$ grid displays the original basis solved from SVD, and the bottom $7\times 2$ grid displays the sparsified basis.
Refer to caption: x5.png

### C.3 Using Differential Invariants in Weak SINDy

In this subsection, we discuss the formulation of weak SINDy and how to implement our strategy of using differential invariants within the weak SINDy framework. To maintain a similar notation to the original works on weak SINDy (Messenger & Bortz, 2021a, b), we use $D_{\alpha_{s}}$ to denote partial derivative operators, where $\alpha_{s}=(s_{1},s_{2},...,s_{p})$ is a multi-index, instead of using subscripts for partial derivatives. Thus, we no longer strictly differentiate subscripts and superscripts–both can be used for indexing lists, vectors, etc.

Given a differential equation in the form

$$ $D_{\alpha_{0}}u=\sum_{s,j}W_{sj}D_{\alpha_{s}}f_{j}(u),$ (50) $$

we can perform integration by parts (i.e., divergence theorem) to move the derivatives from $u$ to some analytic test function and thus bypass the need to estimate derivatives numerically. First, we multiply both sides of ([50](https://arxiv.org/html/2505.12083v1#A3.E50)) by a test function $\phi$ with compact support $\mathcal{B}\subset X$ and integrate over the spacetime domain:

$$ $\int_{X}D_{\alpha_{0}}u(\mathbf{x})\phi(\mathbf{x})d\mathbf{x}=\sum_{s,j}W_{sj }\int_{X}D_{\alpha_{s}}f_{j}(u(\mathbf{x}))\phi(\mathbf{x})d\mathbf{x}$ (51) $$

WLOG, assume that $s_{1}\not=0$, and denote $\alpha_{s^{\prime}}=(s_{1}-1,s_{2},...,s_{p})$. Then, each term in the RHS can be integrated by parts as

$$ $\displaystyle\int_{X}D_{\alpha_{s}}f_{j}(u(\mathbf{x}))\phi(\mathbf{x})d \mathbf{x}$ $\displaystyle=\int_{\mathcal{B}}D_{\alpha_{s}}f_{j}(u(\mathbf{x}))\phi(\mathbf {x})d\mathbf{x}$ (52) $\displaystyle=-\int_{\mathcal{B}}D_{\alpha_{s^{\prime}}}f_{j}(u(\mathbf{x}))D_ {1}\phi(\mathbf{x})d\mathbf{x}+\int_{\partial\mathcal{B}}\nu_{1}D_{\alpha_{s^{ \prime}}}f_{j}(u(\mathbf{x}))\phi(\mathbf{x})d\mathbf{x}$ (53) $\displaystyle=-\int_{\mathcal{B}}D_{\alpha_{s^{\prime}}}f_{j}(u(\mathbf{x}))D_ {1}\phi(\mathbf{x})d\mathbf{x},$ (54) $$

where $D_{1}$ denotes the partial derivative operator w.r.t the first independent variable, and $\nu_{1}$ is the first component of the unit outward normal vector.

Repeating this process until all the derivative operations move from $f_{j}(u)$ to the test function $\phi$, we have

$$ $\int_{X}D_{\alpha_{s}}f_{j}(u(\mathbf{x}))\phi(\mathbf{x})d\mathbf{x}=(-1)^{| \alpha_{s}|}\int_{X}f_{j}(u(\mathbf{x}))D_{\alpha_{s}}\phi(\mathbf{x})d\mathbf {x}$ (55) $$

Similarly for the LHS:

$$ $\int_{X}D_{\alpha_{0}}u(\mathbf{x})\phi(\mathbf{x})d\mathbf{x}=(-1)^{|\alpha_{ 0}|}\int_{X}u(\mathbf{x})D_{\alpha_{0}}\phi(\mathbf{x})d\mathbf{x}$ (56) $$

The final optimization problem is to solve for $\mathbf{b}=\mathbf{G}\mathbf{w}$, where $\mathbf{w}$ is the vectorized coefficient matrix $W$, and each row in $\mathbf{b}$ and $\mathbf{G}$ is given by computing the integrals in ([55](https://arxiv.org/html/2505.12083v1#A3.E55)) and ([56](https://arxiv.org/html/2505.12083v1#A3.E56)) 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 [Section C.2](https://arxiv.org/html/2505.12083v1#A3.SS2), we can enforce symmetry by converting it to a set of linear constraints on 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 $\beta$ 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 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4) about the selected set of symmetry invariants hold. We should note that it is not always possible to find a 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 3.4](https://arxiv.org/html/2505.12083v1#S3.Thmtheorem4). One may ask the following question: can we simply express the equations in terms of differential invariants and apply weak SINDy, similar to [Algorithm 2](https://arxiv.org/html/2505.12083v1#alg2) for the original SINDy formulation? Here, we do not provide a 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 ([50](https://arxiv.org/html/2505.12083v1#A3.E50)), it is more helpful to consider the symmetry invariants as 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 $\{\mathcal{O}_{j}\}$ and a set of fundamental differential invariants $\mathbf{I}=\{\eta_{k}\}$ s.t. every differential invariant can be written as $\mathcal{O}_{j_{1}}...\mathcal{O}_{j_{n}}\eta_{k}$. For the $\mathrm{SO}(2)$ symmetry group in [Example B.5](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem5), one possible choice is

$$ $\eta_{1}=\frac{1}{2}(x^{2}+y^{2}),\ \eta_{2}=u,\ \mathcal{O}_{1}=xD_{y}-yD_{x} ,\ \mathcal{O}_{2}=xD_{x}+yD_{y}.$ (57) $$

We can compose these generating invariant operators to obtain a full library of eligible differential operators up to some order, denoted $\mathbf{D}=\{\mathcal{D}_{j}\}$. The exact compositions can vary and we can choose the most convenient one for subsequent calculations. For the above $\mathrm{SO}(2)$ example, for up to second-order differential operators, we can choose $\{\mathcal{O}_{1},\mathcal{O}_{2},\mathcal{O}_{1}^{2},\mathcal{O}_{2}^{2},
\frac{2}{\eta_{1}}(\mathcal{O}_{1}^{2}+\mathcal{O}_{2}^{2})\}$. Note the last operator is exactly the Laplacian.

Then, the complete set of eligible terms (respecting the symmetry) in the equation is $\{\mathcal{D}_{j}\eta_{k}:\mathcal{D}_{j}\in\mathbf{D},\eta_{k}\in\mathbf{I}\}$. If we assume, as in SINDy, that the governing equation can be written in linear combination of these symmetry invariants, then we can assign a weight for each $\mathcal{D}_{j}\eta_{k}$ and form a coefficient matrix $W=[W_{jk}]$. That is,

$$ $\mathcal{D}_{j_{0}}\eta_{k_{0}}=\sum_{(j,k)\not=(j_{0},k_{0})}W_{jk}\mathcal{D }_{j}\eta_{k}.$ (58) $$

Then, multiplying each side by a test function $\phi(\mathbf{x})$, we have

$$ $\int_{X}\mathcal{D}_{j_{0}}\eta_{k_{0}}\phi(\mathbf{x})d\mathbf{x}=\sum_{(j,k) \not=(j_{0},k_{0})}W_{jk}\int_{X}\mathcal{D}_{j}\eta_{k}\phi(\mathbf{x})d \mathbf{x}.$ (59) $$

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 $D_{\alpha_{s}}$ and ordinary functions $f_{j}(u)$.

To check this, let us explicitly write out the dependency of these operators and fundamental invariants.

#### Case 1

A relatively simple case is when all invariant operators take the form $\mathcal{D}_{j}=\sum_{s}a_{s}(\mathbf{x})D_{\alpha_{s}}$ and $\eta_{k}=\eta_{k}(\mathbf{x},u(\mathbf{x}))$. Each term in the RHS of ([59](https://arxiv.org/html/2505.12083v1#A3.E59)) can be expanded as

$$ $\displaystyle\int_{X}\mathcal{D}_{j}\eta_{k}\phi(\mathbf{x})d\mathbf{x}$ $\displaystyle=\sum_{s}\int_{X}a_{s}(\mathbf{x})D_{\alpha_{s}}\eta_{k}(\mathbf{ x},u(\mathbf{x}))\phi(\mathbf{x})d\mathbf{x}$ (60) $\displaystyle=\sum_{s}(-1)^{|\alpha_{s}|}\int_{X}\eta_{k}(\mathbf{x},u(\mathbf {x}))D_{\alpha_{s}}[a_{s}(\mathbf{x})\phi(\mathbf{x})]d\mathbf{x}$ (61) $$

Evaluating ([61](https://arxiv.org/html/2505.12083v1#A3.E61)) does not require estimating partial derivatives of $u$. Therefore, weak SINDy can be applied to this case quite straightforwardly.

#### Case 2

However, it is not always possible to have all $\mathcal{D}_{j}$ as classical linear differential operators and all $\eta_{k}$ as ordinary functions.
For instance, in [Example B.6](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem6), there are no ordinary symmetry invariants due to the constraint of translation symmetry.

If we still have linear operators $\mathcal{D}_{j}=\sum_{s}a_{s}(\mathbf{x})D_{\alpha_{s}}$, but on the other hand we have differential functions $\eta_{k}=\eta_{k}(\mathbf{x},u^{(n)})$, we can still perform integration by parts as in ([61](https://arxiv.org/html/2505.12083v1#A3.E61)), but the final result becomes

$$ $\sum_{s}(-1)^{|\alpha_{s}|}\int_{X}\eta_{k}(\mathbf{x},u^{(n)})D_{\alpha_{s}}[ a_{s}(\mathbf{x})\phi(\mathbf{x})]d\mathbf{x},$ (62) $$

meaning we still have to evaluate whatever partial derivatives remain in $\eta_{k}$. It is possible that we can 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 3

The most challenging case is when the invariant differential operators explicitly involve the partial derivative, such as $\mathcal{D}_{j}=\sum_{s}a_{s}(\mathbf{x},u^{(n)})D_{\alpha_{s}}$. Then, similar to ([62](https://arxiv.org/html/2505.12083v1#A3.E62)), integration by parts yields:

$$ $\sum_{s}(-1)^{|\alpha_{s}|}\int_{X}\eta_{k}(\mathbf{x},u^{(n)})D_{\alpha_{s}}[ a_{s}(\mathbf{x},u^{(n)})\phi(\mathbf{x})]d\mathbf{x}.$ (63) $$

In this case, we still need to compute the partial derivatives, not only those in $\eta_{k}$, but also those arising from $a_{s}$ and $D_{\alpha_{s}}(a_{s})$. The latter might involve higher-order derivatives and the benefit of using the weak formulation may further diminish.

## Appendix D Additional Experiment Results

Contents of this section include:

- •
[Section D.1](https://arxiv.org/html/2505.12083v1#A4.SS1): Results for some variants of the sparse regression models considered in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1).
- •
[Section D.2](https://arxiv.org/html/2505.12083v1#A4.SS2): Results for genetic programming-based algorithms under different computational budgets.
- •
[Section D.3](https://arxiv.org/html/2505.12083v1#A4.SS3): Samples of equations discovered by different methods.
- •
[Section D.4](https://arxiv.org/html/2505.12083v1#A4.SS4): Visualized prediction errors of equations discovered by different methods.

### D.1 Variant Sparse Regression Models

**Table 3: Results of sparse regression models on the Boussinesq equation and the reaction-diffusion system. $\mathcal{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 1](https://arxiv.org/html/2505.12083v1#S4.T1).**
| Method | Boussinesq ([7](https://arxiv.org/html/2505.12083v1#S4.E7)) | Reaction-diffusion ([10](https://arxiv.org/html/2505.12083v1#S4.E10)) |  |  |
| --- | --- | --- | --- | --- |
| $\mathcal{C}$ $\downarrow$ | SP $\uparrow$ | $\mathcal{C}$ $\downarrow$ | SP $\uparrow$ |  |
| PySINDy | 15 | 0.00 | 38 | 0.53 |
| PySINDy^∗ | 21 | 1.00 | 468 | 0.00 |
| SI | 13 | 1.00 | 28 | 0.54 |
| SI-aligned | - | - | 14 | 0.56 |

The original implementation of PySINDy (de Silva et al., 2020; Kaptanoglu et al., 2022) does not allow functions to be applied to partial derivative terms. As a result, terms such as $u_{x}^{2}$ cannot be modeled. This leads to its failure to discover the Boussinesq equation ([7](https://arxiv.org/html/2505.12083v1#S4.E7)), as we have shown in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1).

We modify the implementation and include an additional set of results with different libraries, denoted as PySINDy^∗ in [Table 3](https://arxiv.org/html/2505.12083v1#A4.T3).
The PySINDy^∗ model supports a wider range of library functions, including functions of partial derivatives, e.g., $u_{x}^{2}$. A complete description of the hypothesis spaces of different sparse regression-based methods is available in [Section E.5](https://arxiv.org/html/2505.12083v1#A5.SS5).

As [Table 3](https://arxiv.org/html/2505.12083v1#A4.T3) shows, PySINDy^∗ succeeds in the Boussinesq equation. However, it fails in the reaction-diffusion system because its parameter space becomes too large due to a higher-dimensional total space $X\times U\simeq\mathbb{R}^{2}\times\mathbb{R}^{2}$. This augments the point that SINDy’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 is sure to fail. If the library is too large, the optimization problem becomes more difficult in the high-dimensional parameter space.
On the other hand, by introducing the inductive bias of symmetry, our method automatically identifies a proper function library that contains all the necessary terms for a PDE with a specific symmetry group, but not other redundant terms.

We include another model in [Table 3](https://arxiv.org/html/2505.12083v1#A4.T3), SI-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 in a constrained subspace) as in 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 3.3](https://arxiv.org/html/2505.12083v1#S3.SS3) and [Section C.2](https://arxiv.org/html/2505.12083v1#A3.SS2). 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 3](https://arxiv.org/html/2505.12083v1#A4.T3) for reference.

For the reaction-diffusion system, SI-aligned has a $14$-dimensional parameter space. The basis for its parameter space is visualized in [Figure 4](https://arxiv.org/html/2505.12083v1#A3.F4). It achieves a slightly higher success probability than SI (regression with symmetry invariants) and PySINDy (without symmetry 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 SI (regression with symmetry invariants) by multiplying both sides by $u_{x}^{2}$.

We note that the results on the reaction-diffusion system in [Table 3](https://arxiv.org/html/2505.12083v1#A4.T3) are for models with the original SINDy formulation, in contrast to the weak SINDy formulation used in [Figure 3](https://arxiv.org/html/2505.12083v1#S4.F3). Therefore, the results in [Figure 3](https://arxiv.org/html/2505.12083v1#S4.F3) should not be directly compared to those in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1) and [Table 3](https://arxiv.org/html/2505.12083v1#A4.T3).

### D.2 Genetic Programming

Figure: Figure 5: Success Probabilities of GP-based methods on different systems. Our method with symmetry invariants can discover the correct equations with fewer iterations.

For each system in [Section 4.1](https://arxiv.org/html/2505.12083v1#S4.SS1), we run the genetic programming discovery algorithm with three different iteration counts, but otherwise keep all hyperparameters constant. In Figure [5](https://arxiv.org/html/2505.12083v1#A4.F5), we plot the success probability as a function of the iteration count for both the base GP algorithm and our method that uses symmetry invariants.

In all cases, we find that using symmetry invariants results in a higher success probability in comparison to unmodified PySR. Specifically, for the Boussinesq equation, our method achieves a 100% chance of discovery with 5 iterations, whereas even with 3 times the number of iterations, PySR only yields a 90% success probability. This highlights that using invariants improves the efficiency of equation discovery. For Darcy flow and Reaction-Diffusion, we find that the base genetic programming algorithm fails to ever make a correct prediction. On the other hand, using symmetry invariants leads to a successful discovery the majority of the time.

We finally note that increasing the number of iterations to $200$ for Darcy flow slightly lowers the success probability when using symmetry invariants. We hypothesize this is because at higher iterations, the search process begins to overfit and introduces extraneous low-order terms. While we already drop some terms with small enough coefficients, future works may consider a more refined filtration process.

### D.3 Samples of Discovered Equations

In [Table 4](https://arxiv.org/html/2505.12083v1#A4.T4), we list some randomly selected equations discovered by different methods for the Boussinesq equation ([7](https://arxiv.org/html/2505.12083v1#S4.E7)). Some methods almost consistently discover correct/incorrect equations (i.e., have success probabilities close to 1 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: a correct equation and an incorrect one.

The ground truth equation in the original variables is given in ([7](https://arxiv.org/html/2505.12083v1#S4.E7)). The ground truth equation in the symmetry invariants is given by

$$ $\eta_{(0,2)}+\eta_{(0,0)}\eta_{(2,0)}+\eta_{(4,0)}+1=0$ (64) $$

**Table 4: Samples of discovered equations from the observed solution of the Boussinesq equation ([7](https://arxiv.org/html/2505.12083v1#S4.E7)). For GP-based methods, we include results from different numbers of iterations (indicated by "$N$ its"). For transformer-based methods, we include two samples for each method because of the large variance of discovered equations from different runs.**
| Method | Equation sample(s) |  |
| --- | --- | --- |
| Sparse regression | PySINDy | $u_{tt}=-1.01u_{xxxx}-0.79uu_{xx}$ |
| PySINDy^∗ | $u_{tt}=-1.01u_{xxxx}-0.99u_{x}^{2}-0.98uu_{xx}$ |  |
| SI | $\eta_{(0,2)}=-1.00-1.00\eta_{(4,0)}-1.00\eta_{(0,0)}\eta_{(2,0)}$ |  |
| Genetic programming | PySR (5 its) | $uu_{xx}+1.00u_{tt}+u_{xxxx}=0$ |
| PySR (15 its) | $uu_{xx}+u_{tt}+u_{x}^{2}+1.00u_{xxxx}=0$ |  |
| SI (5 its) | $1.00\eta_{(0,0)}\eta_{(2,0)}+1.00\eta_{(0,2)}+1.00\eta_{(4,0)}+1=0$ |  |
| Transformer | E2E | (1) $u_{tt}=-1.13uu_{xx}-0.98u_{xxxx}-0.30|u_{x}|$ |
| (2) $u_{tt}=-0.85uu_{xx}-0.75u_{x}^{2}-0.99u_{xxxx}$ |  |  |
| SI | (1) $\eta_{(0,2)}=-1.05\eta_{(0,0)}\eta_{(2,0)}-1.00\eta_{(4,0)}-0.96$ |  |
| (2) $\eta_{(0,2)}=-0.81\eta_{(0,0)}\eta_{(2,0)}-0.40\eta_{(0,0)}-0.98\eta_{(4,0)}-0<br>.90$ |  |  |

[Table 5](https://arxiv.org/html/2505.12083v1#A4.T5) lists the equation samples discovered from the Darcy flow dataset. The ground truth equation in original variables is given in ([8](https://arxiv.org/html/2505.12083v1#S4.E8)), and the ground truth equation in symmetry invariants is given by

$$ $8\zeta_{2}-\Delta u-e^{4R^{2}}=0,$ (65) $$

where $\zeta_{2}=xu_{x}+yu_{y}$, $\Delta u=u_{xx}+u_{yy}$, and $R^{2}=x^{2}+y^{2}$ are among the rotational invariants used in symbolic regression.

**Table 5: Samples of discovered equations for the Darcy flow dataset.**
| Method | Equation sample |  |
| --- | --- | --- |
| Genetic programming | PySR | $u-0.47x^{2}y^{2}-0.38e^{0.09(u_{xx}+u_{yy})}+0.20=0$ |
| SI | $\zeta_{2}-0.13\Delta u-0.13e^{4.01R^{2}}=0$ |  |
| Transformer | E2E | $u_{xx}=-7.43\sqrt{u^{2}+0.65u_{x}^{2}}$ |
| SI | $\Delta u=-2.56u+0.85\zeta_{2}+0.29$ |  |

Finally, [Table 6](https://arxiv.org/html/2505.12083v1#A4.T6) lists the equation samples discovered from the reaction-diffusion dataset. The ground truth equation in original variables is given in ([10](https://arxiv.org/html/2505.12083v1#S4.E10)) with $d_{1}=d_{2}=0.1$, and the ground truth equation in symmetry invariants is given by

$$ $\displaystyle I_{t}$ $\displaystyle=0.1(I_{xx}+I_{yy})+A(1-A)$ (66) $\displaystyle E_{t}$ $\displaystyle=0.1(E_{xx}+E_{yy})-A^{2}$ (67) $$

where $I_{\mu}=uu_{\mu}+vv_{\mu}$ and $E_{\mu}=-vu_{\mu}+uv_{\mu}$ for any multiindex $\mu$ of $t,x,y$, and $A=u^{2}+v^{2}$.

**Table 6: Samples of discovered equations for the reaction-diffusion system dataset. Each discovered result contains two equations, since this is an evolution system with two dependent variables $u,v$.**
| Method | Equation sample |  |
| --- | --- | --- |
| Sparse regression | PySINDy | $\begin{cases}u_{t}=0.96u-0.97u^{3}+1.00^{3}-0.97uv^{2}+1.00u^{2}v+0.09u_{xx}+0<br>.09u_{yy}\\<br>v_{t}=0.96v-1.00u^{3}-0.97v^{3}-1.00uv^{2}-0.96u^{2}v+0.09v_{xx}+0.09v_{yy}<br>\end{cases}$ |
| PySINDy^∗ | $\begin{cases}u_{t}=0.21u-0.24u^{3}+1.00v^{3}-0.23uv^{2}+0.99u^{2}v\\<br>v_{t}=0.21v-1.01u^{3}-0.24v^{3}-0.99uv^{2}-0.23u^{2}v\\<br>\end{cases}$ |  |
| SI | $\begin{cases}I_{t}=0.10I_{xx}+0.10I_{yy}+0.96A-0.96A^{2}\\<br>E_{t}=0.10E_{xx}+0.10E_{yy}-1.00A^{2}\end{cases}$ |  |
| SI-aligned | $\begin{cases}u_{t}=0.95u-0.96u^{3}+1.00v^{3}-0.96uv^{2}+1.00u^{2}v+0.09u_{xx}+<br>0.09u_{yy}\\<br>v_{t}=0.95v-1.00u^{3}-0.96v^{3}-1.00uv^{2}-0.96u^{2}v+0.09v_{xx}+0.09v_{yy}<br>\end{cases}$ |  |
| Genetic programming | PySR | $\begin{cases}u_{t}=0.92v\\<br>v_{t}=-0.92u\\<br>\end{cases}$ |
| SI | $\begin{cases}I_{t}=0.10I_{xx}+0.10I_{yy}+A-1.00A^{2}\\<br>E_{t}=0.10E_{xx}+0.10E_{yy}-1.00A^{2}\end{cases}$ |  |
| Transformer | E2E | $\begin{cases}u_{t}=0.89u_{y}\\<br>v_{t}=-0.91u\end{cases}$ |
| SI | $\begin{cases}I_{t}=0\\<br>E_{t}=0.50\arctan(0.45E_{y}-0.31E_{y}/(-540.12AE_{y}+...)+...)+...\end{cases}$ |  |

### D.4 Prediction Errors of Discovered Equations

In [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1), 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, e.g., $t\in[0,20]$ for the Boussinesq equation and $t\in[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.

Figure: Figure 6: Prediction error over time using the discovered equations.
Refer to caption: x10.png

In addition to the prediction error at the end of the simulation time, [Figure 6](https://arxiv.org/html/2505.12083v1#A4.F6) 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 6](https://arxiv.org/html/2505.12083v1#A4.F6) are consistent with those in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1). Generally, the discovered equations with smaller prediction errors at the end of the simulation time also have lower prediction errors throughout the entire time interval.

For Darcy flow ([8](https://arxiv.org/html/2505.12083v1#S4.E8)), since it describes the steady state of a system and does not involve time derivatives, we do not simulate the discovered PDEs. Instead, we evaluate each discovered PDE $F(\mathbf{x},u^{(n)})=0$ on the test dataset $\{(\mathbf{x},u^{(n)}):\mathbf{x}\in\Omega\}$ and report the residual as the prediction error. In addition to the average error over all the spatial grid points reported in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1), we visualize the error heatmaps over the grid in [Figure 7](https://arxiv.org/html/2505.12083v1#A4.F7). It can be observed that the discovered equations with symmetry invariants have lower errors across the entire grid.

Figure: Figure 7: Prediction error of discovered equations from genetic programming methods for Darcy flow. Left: genetic programming with regular variables. Right: genetic programming with symmetry invariants.
Refer to caption: x12.png

## Appendix E Experiment Details

In this section, we describe the experiment setups required to reproduce the experiments. In terms of computational resources, our experiments are conducted with 12 INTEL(R) XEON(R) PLATINUM 8558 CPUs and should be reproducible within minutes with any modern CPUs.

### E.1 Data generation

#### Boussinesq equation

The equation is solved using a Fourier pseudospectral method for spatial derivatives and a fourth-order Runge-Kutta (RK4) scheme for time integration. The solution is computed on a periodic spatial domain $[-L,L]$ with $N=256$ grid points. The equation is reformulated as a first-order system in time by introducing $v=u_{t}$, 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 $u$ are recorded in the dataset for equation discovery. The simulation starts from an initial condition of $u(x)=0.5e^{-x^{2}}$ and $u_{t}=0$ and proceeds up to a final time $T=20$ with a time step of $\Delta t=0.001$. Starting from the solution at $T=20$, we simulate for another $T^{\prime}=20$ with the same configuration to obtain a test dataset for evaluating prediction errors of the discovered equations.

#### Darcy flow

We use the data generation code(^2^22[https://github.com/pdebench/PDEBench/tree/main/pdebench/data_gen/data_gen_NLE/ReactionDiffusionEq](https://github.com/pdebench/PDEBench/tree/main/pdebench/data_gen/data_gen_NLE/ReactionDiffusionEq))
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

$$ $u_{t}(\mathbf{x},t)-\nabla(a(\mathbf{x})\nabla u(\mathbf{x},t))=f(x),\mathbf{x }\in(-0.5,0.5)^{2},$ (68) $$

with $a(\mathbf{x})=e^{-4\|\mathbf{x}\|^{2}_{2}}$ and $f(\mathbf{x})=1$.

#### Reaction-diffusion

We use the data generation code(^3^33[https://github.com/dynamicslab/pysindy/blob/master/examples/10_PDEFIND_examples.ipynb](https://github.com/dynamicslab/pysindy/blob/master/examples/10_PDEFIND_examples.ipynb)) from PySINDy (de Silva et al., 2020; Kaptanoglu et al., 2022). The spatial domain is $[-10,10]\times[-10,10]$ with $128$ grid points in each direction. The simulation proceeds up to a final time $T=10$ with a time step $\Delta t=0.05$. We perturb the numerical solution by a $0.05\%$ noise and record the values of $u,v$ to the dataset for equation discovery. Starting from the solution at $T=10$, we simulate for another $T^{\prime}=10$ with the same configuration to obtain a test dataset for evaluating prediction errors of the discovered equations.

### E.2 Sparse regression

#### Boussinesq equation

For SINDy with original variables, we fix $u_{tt}$ as the LHS of the equation and include functions of up to $4$th-order derivatives on the RHS. For PySINDy in [Table 1](https://arxiv.org/html/2505.12083v1#S4.T1), the library contains monomials on $U^{(4)}$ with degree in $u$ no larger than $2$ and degree in any partial derivative terms $u_{\alpha}$ no larger than $1$. For example, $u^{2}u_{x}$ is included, but $u^{3}$, $u_{x}^{2}$ are not. For PySINDy^∗, the library contains all monomials on $U^{(4)}$ up to degree $2$. For example, $u_{x}^{2}$ and $uu_{x}$ are included. Note that the PySINDy^∗ library does not contain all functions in the original PySINDy library, e.g., $u^{2}u_{x}$ is not included because it has degree $3$.

Our method, SI, uses the invariant set in [Example B.6](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem6) for sparse regression. Specifically, $\eta_{(0,2)}=u_{tt}/u_{x}^{2}$ 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_{x}$ on the denominator, we remove the data points with small $|u_{x}|$ to avoid numerical issues.

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 SI-raw is performed after subsampling. The threshold value for sequential thresholding is set to $0.25$, and the coefficient for $L_{2}$ regularization is set to $0.05$.

#### Darcy flow

Sparse regression-based methods are not directly applicable to Darcy flow ([8](https://arxiv.org/html/2505.12083v1#S4.E8)) because there exist terms such as $e^{-4(x^{2}+y^{2})}$. While it is still possible to include all necessary terms in the function library so that the equation can be written in the linear combination form ([4](https://arxiv.org/html/2505.12083v1#S2.E4)), 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_{t}$ and $v_{t}$ as 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 $u,v$ up to degree $3$ and all spatial derivatives up 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 $\{t,x,y,u^{2}+v^{2}\}\bigcup\{\mathbf{u}\cdot\mathbf{u}_{\mu}\}\bigcup\{
\mathbf{u}^{\perp}\cdot\mathbf{u}_{\mu}\}$, where $\mathbf{u}=(u,v)^{T}$ and $\mu$ is a multiindex of $t,x,y$. We will denote $I_{\mu}=\mathbf{u}\cdot\mathbf{u}_{\mu}$ and $E_{\mu}=\mathbf{u}^{\perp}\cdot\mathbf{u}_{\mu}$. We use $I_{t}$ and $E_{t}$ 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 $L_{2}$ regularization is set to 0 0 for SINDy with original variables and $0.1$ for our method with symmetry invariants.

For the experiments with different levels of noise ([Section 4.4](https://arxiv.org/html/2505.12083v1#S4.SS4)), we use weak SINDy as the base algorithm. The function library is the same as PySINDy as described above. To enforce symmetry, instead of directly using the symmetry invariants, we derive a set of linear constraints on the sparse regression parameters to adapt to weak SINDy. This procedure is further described in [Section C.3](https://arxiv.org/html/2505.12083v1#A3.SS3).

### E.3 Genetic Programming

In all experiments, to determine if an equation matches the ground truth we first expand the prediction into a sum of monomial terms. We then eliminate all terms whose relative coefficient is below $0.01$. For each 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 a single equation, PySR finally produces a hall-of-fame that consists of multiple candidate solutions with varying complexities. To finally pick a single prediction, we use a selection strategy equivalent to the “best” option from PySR.

#### Boussinesq equation

For the Boussinesq equation ([7](https://arxiv.org/html/2505.12083v1#S4.E7)), 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 each variable in $(\mathbf{x},u^{(4)})$ 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 B.6](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem6) and the RHS as a 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 Darcy experiment ([8](https://arxiv.org/html/2505.12083v1#S4.E8)), 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 $(\mathbf{x},u^{(2)})$ for the LHS and the RHS is then a function of the unused variables. Likewise when using invariants, we fix the LHS to each possible invariant specified in [Example B.5](https://arxiv.org/html/2505.12083v1#A2.Thmtheorem5) 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 ([10](https://arxiv.org/html/2505.12083v1#S4.E10)), 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, to have $127$ populations of size $64$, and to have a fraction-replaced coefficient of $0.5$.

In the ordinary variable case, we fix the LHS as either $u_{tt}$ or $v_{tt}$ and allow the RHS to be a function of all other variables in $(\mathbf{x},u^{(2)})$. When using invariants, the LHS is fixed to be either $I_{t}$ or $E_{t}$ 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(^4^44[https://github.com/facebookresearch/symbolicregression/blob/main/Example.ipynb](https://github.com/facebookresearch/symbolicregression/blob/main/Example.ipynb)) from Kamienny et al. (2022). The transformer-based 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.

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 $u_{tt}$ for original variables and $\eta_{(0,2)}$ for symmetry invariants. For the Darcy flow, the LHS is fixed to $u_{xx}$ for original variables and $\Delta u$ for symmetry invariants. For the reaction-diffusion system, the LHS is fixed to $u_{t},v_{t}$ for original variables and $I_{t},E_{t}$ for symmetry invariants.

### E.5 Hypothesis Spaces of Equation Discovery Algorithms

[Table 7](https://arxiv.org/html/2505.12083v1#A5.T7) and [Table 8](https://arxiv.org/html/2505.12083v1#A5.T8) describe the hypothesis spaces of different equation discovery algorithms when applied to the Boussinesq equation and the reaction-diffusion system.

**Table 7: Hypothesis spaces of different equation discovery algorithms for the Boussinesq equation.**
| Method | Hypothesis space |  |
| --- | --- | --- |
| Sparse Regression | PySINDy | $u_{tt}=W\bm{\theta}(u^{(4)}),\{\theta^{j}\}=\{ab:a\in\mathrm{Mono}_{\leq 2}(U)<br>,b\in\{1,u_{x},...,u_{xxxx}\}\}$ |
| PySINDy^∗ | $u_{tt}=P(u^{(4)})\in\mathrm{Poly}_{\leq 2}(U^{(4)})$ |  |
| SI | $\eta_{(0,2)}=P(\bm{\eta})\in\mathrm{Poly}_{\leq 2}(\{\eta_{(\alpha,\beta)}\}<br>\backslash\{\eta_{(0,2)}\})$ |  |
| Genetic Programming | PySR | $z^{j}=f(\mathbf{z}^{-j})$ for $\mathbf{z}=(\mathbf{x},u^{(4)})$ and some $j$ |
| SI | $\eta_{(\alpha_{0},\beta_{0})}=f(\bm{\eta}_{-(\alpha_{0},\beta_{0})})$ for $\bm{\eta}=\{\eta_{(\alpha,\beta)}:\alpha+\beta\leq 4\}$ and some $(\alpha_{0},\beta_{0})$ |  |

**Table 8: Hypothesis spaces of different equation discovery algorithms for 2D reaction-diffusion. $\mathbf{u}^{(n)}\in U^{(n)}$ denotes the collection of all up to $n$th order spatial derivatives. $\alpha=[\alpha_{1},\alpha_{2}]$ is the multiindex for spatial variables. $\mathbf{x}=(x,y,t)$. $A=u^{2}+v^{2}$.**
| Method | Hypothesis space |  |
| --- | --- | --- |
| Sparse Regression | PySINDy | $\mathbf{u}_{t}=W\bm{\theta}(\mathbf{u}^{(2)}),\{\theta^{j}\}=\mathrm{Mono}_{<br>\leq 3}(U)\bigcup\{\mathbf{u}_{\alpha}:|\alpha|\leq 2\}$ |
| PySINDy^∗ | $\mathbf{u}_{t}=P(\mathbf{u}^{(2)})\in\mathrm{Poly}_{\leq 3}(U^{(2)})$ |  |
| SI | $[I_{t},E_{t}]^{T}=P\in\mathrm{Poly}_{\leq 2}(A,\mathbf{x},I_{\alpha},E_{\alpha<br>};|\alpha|\leq 2)$ |  |
| SI-aligned | $\mathbf{u}_{t}=W\bm{\theta}(\mathbf{u}^{(2)}),W^{jk}=Q^{ijk}\beta^{i}$ for some precomputed $Q$ |  |
| Genetic Programming | PySR | $\mathbf{u}_{t}=\mathbf{f}(\mathbf{x},\mathbf{u}^{(2)})$ |
| SI | $[I_{t},E_{t}]^{T}=\mathbf{f}(A,\mathbf{x},I_{\alpha},E_{\alpha};|\alpha|\leq 2)$ |  |

## Appendix F Broader Impacts

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.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 4
presentation: 4
questions: - Could you quantify the human effort involved in deriving and "adjusting" the invariants for a system? What about when the invariants become so complex to be uninterpretable, thus negating the benefit of symbolic regression?
- Algorithm 1 proposes iterating through each invariant to serve as LHS and selecting the one that minimizes regression error. This seems computationally intensive and potentially non-physical. Could this approach not lead to discovering a mathematically valid but physically meaningless rearrangement of the true equation (like expressing $u$ as a complex function of $u_t$ and $u_xx$)?
- For the imperfect symmetry experiments, how is the regularization strength on the symmetry-breaking terms chosen ($d_1$ and $d_2$)? These are sensitive hyperparameters that would require careful tuning on a validation set, potentially limiting the method's utility when ground truth is unavailable.
- How is the "specific timestep" of the computation of the prediction errors (PE) chosen?
- Could you elaborate on the choice of median RMSE instead of averages due to "tremendous prediction errors"? Are these outliers vaues, or happen more frequently? Histograms of these estimates could further convey the effectiveness of the method, since it has a success probability >54% (when applicable). 
- Concerning GP and Appx. D2, why did you set up to 15 iterations for Boussinesq, and hundreds for the other cases? Considering that the default is 100 and performance was already on par with SI at 15. Do you have any intuition why PySR fails completely on Darcy and RD?
- Is there a reason why you picked RD and SINDy for the experiments in 4.4, while not investigating the Boussinesq case with noisy data and other baselines? As you also mentioned, SINDy is not the best approach when dealing noisy data. Additionally, in the Boussinesq experiment, the baseline PySINDy fails because its default library lacks the necessary terms, but a modified PySINDy* with said terms succeeds perfectly (Table 3 in the appendix). This suggests the baseline's failure is one of library design, not a fundamental flaw. The key advantage of the invariant method is that it "automates" the creation of a correct library, and this point should be made more transparent, rather than presenting a 0% vs. 100% success rate against a limited baseline (Table 1 in the main text).
soundness: 3
strengths: This paper is well-writen and mathematically elegant. The main strengths are:
- Its primary contribution, of integrating symmetry via differential invariants, is a sound and rigorous approach based in established mathematical theory.
- This framework is presented as an "enhancement" that can be applied to the main SR methods, making this contribution potentially more impactful.
- Within the context of the chosen experimental settings, the methods show remarkable effectiveness compared to the baselines. There are also studies that address noisy conditions and imperfect symmetries, amenable for real-world cases.
- Clear Figures, Proofs, reproducible code, and a broad, extended appendix further improve the quality and soundness of this scientific work as a whole.
summary: This work proposes a framework for incorporating known physical symmetries in the process of PDE equation discovery via symbolic regression (SR) from data. The central idea is to leverage the mathematical theory of Lie groups and differential invariants to enhance the symbolic regression algorithms. These invariants are functions of the original variables (and their derivatives) that by definition remain unchanged under the symmetry transformations, hence the resulting discovered equation is guaranteed to respect these symmetries by construction, and are more easily discovered. The authors further test this approach by integrating it with three of the main classes of SR (SINDy, PySR, and symbolic transformers), showing improvements in success rates and predictions on three canonical PDE systems.

Although the paper has a novel approach to enforce symmetries in data for the challenging task of SR, is very well written and clear, with an extensive appendix and well-documetned code, its novelty and experimental evaluations are limited. Most of the presented theory is based on Olver (1993), and in practice their method is effective on a few specific cases (ref. Tab 1). Most notably, this approach requires that the symmetry group is already given and differential invariants have to be computed and adapted to the SR methods (when feasible). These limitations are discussed briefly only at the end of Sec. 5, but I believe this should be clearly stated earlier in the presentation.
This work deserves attention from the ML community, and while this venue may not be the most appropriate (even due to its page limit), this line of research is worth more investigation and discussion. Also, the contribution on the computation of high-order differential invariant (Prop 3.3) may be worth a standalone publication in a mathematical venue or journal.
Finally, I am willing to discuss with the authors, AC and other reviewers, and reconsider my evaluation.
weaknesses: Despite the strengths, the paper suffers from significant limitations in its core assumptions, novelty, and experimental scope, which are not given sufficient weight in the main body of the text.

-  The most critical weakness is the core assumption that the symmetry group of the PDE is known a priori. The paper positions itself as a tool for equation discovery, but, as stated, in many scientific discovery scenarios the underlying symmetries are unknown and are part of what needs to be discovered. Besides all the benchmarks and comparisons that are necessary, it would have been more interesting to address this issue and limitation. Experimental tests on black-box, symbolic-unknown real-world dataset on which some symmetries are assumed would highly showcase the effectiveness of this approach, and broaden its impact.
- While the application of this approach to a broad set of SR methods is a solid contribution to the ML community, the core conceptual idea of using symmetries is a standard approach in theoretical physics. The math presented in Sec 3.1 also belongs to the Background and is not novel; the noteworthy mathematical contributions are Prop 3.3 and 3.4, and in practice there are some limitations (could you elaborate on line 248 "We evaluate invariants on the dataset only where they are well-defined"?). The paper's claim to be "the first attempt to strictly enforce general symmetries [...] for general symbolic regression methods" (lines 759) hinges heavily on the word "general." The conceptual leap seems less significant than presented, and the work can be viewed as an incremental, yet valuable, extension of prior ideas into a broader software framework.
- The paper significantly downplays the practical difficulty of deriving the differential invariants. The authors state that solving for invariants "is easy with any symbolic computation package" but later concede that the results "may be complicated and require ad-hoc adjustment for better interpretability and compatibility" (lines 481-484). This "one-time effort" can be extremely non-trivial for complex symmetries or high-dimensional systems, requiring substantial domain expertise. This represents a major "pre-processing" step that acts as a significant barrier to entry for practitioners, and undermines the goal of automated equation discovery.
- Lastly, the experimental setting, while clean, is limited in scope and does not convincingly demonstrate robustness for real-world challenges or more complex systems.

Minor comments that did not affect the score:
- The compared baseline methods, although established in the literature, have known limitations, and other more powerful algorithms are available. See SR Bench (https://cavalab.org/srbench/) or TPSR (Shojaee 2023) instead of E2E. 
- Error estimates on your results are missing, making the comparison of the results less valuable. You could use a simple bootstrap of median values as an approximation of median values, or report histograms of the outcomes.
- For clarity, Appendix D and E should be reversed, as they appear in that order in the text.
- There are some arXiv references, which should be avoided as they are not peer-reviewed. Refer to peer-reviewed venues instead (if applicable).

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. For higher $p,n$, the count of derivatives explodes. What is the computational complexity (symbolic time) of producing invariants via $\text{pr}^{(n)}v(\eta)=0$ and the Prop. 3.3 recursion? Any timings and symmetry invariant size for your three PDEs in Table 1?

2. Algorithm 1 tries each invariant as LHS and picks the lowest relative error. How do you mitigate scale effects or overfitting to noise in this selection?
soundness: 2
strengths: 1. **Significance**: The paper focuses on compressing the vast and high-dimensional search space for PDEs, a central problem in existing PDE discovery methods.

2. **Clarity:** The paper has rigorous and well-defined notations and terminologies, with a concrete definition of the construction of differential invariants, algorithm 1 for explicit SR with invariants, and appendix samples/tables making the pipeline replicable and interpretable.

3. **Adaptive Design to Several Baselines:** Adapt symmetry invariants into the major baselines (SINDy, GP, and Transformer), making the proposed method flexible to different problem settings, and showing a consistent improvement for accuracy and the complexity management capability among those baselines.
summary: This paper introduces a general framework to incorporate the differential invariants of partial differential equations (PDEs) derived from Lie symmetry theory into several prevailing symbolic regression methodologies. Experiments on several PDE systems and noise settings verify the effectiveness of the proposed symmetry invariants for more accurate, robust, and symmetry-consistent PDE discovery.
weaknesses: 1. **Lack of Statistical Significance**: The paper claims the Table 1 results are based on 100 runs for each algorithm; however, only the median of PE is reported. A confidential interval can strengthen the statistical significance of the results.

2. **Lack of Computational Efficiency Analysis**:  The paper claims in line 21 that the proposed approach can "improve their (sparse regression's and genetic programming's) accuracy and *efficiency*." However, only accuracies (SP and PE) are reported in Table 1 and Figure 3. There is no efficiency analysis, such as running time or search space size benchmarking for the experiments in the main context, to support the claim.

---

confidence: 4
contribution: 1
initial_rating: 2
presentation: 2
questions: - For the Reaction-Diffusion equation, why does the method consider the SO(2) group acting on the dependent variables (u, v) rather than on the independent variables (x, y)? Where does this prior knowledge originate? In fact, for this system, the SO(2) symmetry on (x, y) is more general, while the SO(2) symmetry on (u, v) only holds when $d_1=d_2$—it is a case-specific condition.

- If the assumed symmetry is incorrect, to what extent does the proposed method fail? This will be an important issue when integrating the approach with symmetry discovery in future work.

- In the case of the 2D Navier–Stokes equations, how does the method handle the continuity constraint $\nabla \cdot u = 0$? How is the number of equations determined?

- How can discrete symmetries be incorporated into the proposed framework?
soundness: 2
strengths: - From the perspective of differential invariants, this paper provides an efficient method for discovering PDEs when symmetries are known.

- The proposed method does not alter the existing symbolic regression framework except for the input variables, making it easy to implement.

- The experimental setting is reasonably comprehensive, including evaluations under noisy data and imperfect symmetry scenarios.
summary: This paper proposes a method for discovering symbolic differential equations from data using symmetry invariants.
The authors argue that by employing the differential invariants of symmetry groups as fundamental variables in symbolic regression, the discovered equations can be guaranteed to satisfy specified symmetries, thereby narrowing the search space and improving discovery efficiency and accuracy.
The method is integrated into existing symbolic regression algorithms such as sparse regression (e.g., SINDy) and genetic programming, and is validated on several PDE systems (e.g., Boussinesq equation, Darcy flow, reaction-diffusion systems), demonstrating its robustness in the presence of noisy data and imperfect symmetries.
weaknesses: - Clarification of contributions. The methodology section of this paper is organized into four main components: the definition of symmetry invariants (Section 3.1), the computation of symmetry invariants (Section 3.2), the preprocessing of invariants as variables for the symbolic regression algorithm (Section 3.3), and the constraint relaxation for systems with imperfect symmetry (Section 3.4).
The framework largely builds upon existing techniques rather than introducing a substantially new theoretical formulation.
In particular, the approach for computing higher-order differential invariants is already well-established (e.g., https://doi.org/10.1007/978-1-4684-0274-2 and https://arxiv.org/abs/2307.05432v2);
yet the authors dedicate an entire section (Section 3.2) to this topic without appropriate citations or clarification of novelty, which makes it difficult for readers to identify the paper’s specific contributions.
In addition, recent studies (e.g., https://arxiv.org/abs/2402.03747
) have also explored the integration of symmetry and PDE discovery.
A more comprehensive discussion and comparison with these related works would help clarify how this paper advances beyond prior efforts.

- Non-uniqueness of fundamental invariants.
Even for the same symmetry group, a complete set of fundamental (ordinary) invariants is generally not unique.
Different derivation procedures or generator choices can lead to distinct but equivalent invariant forms.
In Appendix B.4, the authors directly state that the two ordinary invariants for SO(2) are $\tfrac{1}{2}(x^2 + y^2)$ and $u$ without derivation.
While this choice is straightforward for simple symmetry groups, it becomes ambiguous for more complex transformations (see Chapter 2.4 of https://doi.org/10.1007/978-1-4684-0274-2 for more examples).
**If the selected invariants do not align with those appearing in the ground-truth PDE, the symbolic regression process may fail to recover the correct equation.**
Furthermore, differences in the derivation of ordinary invariants can propagate to the higher-order invariants computed via Proposition 3.3, which serve as meta-variables in SR methods, thereby introducing instability and inconsistency in more complex scenarios.

- The experimental section of the paper only considers common and relatively simple symmetry groups such as SO(2), scaling, and translation, which does not sufficiently validate the effectiveness of the proposed method.
More importantly, the approach relies on a strong assumption—that explicit expressions for the infinitesimal generators of the symmetry group are known in advance.
However, this requirement is difficult to satisfy in most real-world scenarios, as the Lie point symmetries of PDEs can have highly complex and nonlinear forms.
The authors mention in Section 5 that future work could integrate symmetry discovery methods, but most existing approaches involve neural network modules that are not analytically interpretable, making it impractical to obtain explicit generator expressions (e.g., the encoder–decoder in LaLiGAN (Yang et al., 2024) or the MLP used by Ko et al., 2024).
Consequently, it remains unclear whether the proposed method can genuinely facilitate the discovery of entirely unknown PDEs from data, rather than merely rediscovering toy examples already well understood by human experts.

- None of the experimental results include error bars or measures of variability, making it difficult to fully assess the reliability, robustness, and fairness of the reported findings.

- Inconsistency between theory, algorithm, and experiments.
Proposition 3.3 assumes the presence of only one prolonged vector field, which is quite limited and does not hold for more complex symmetry groups such as SO(3) or the Lorentz group that involve multiple generators.
In contrast, Algorithm 1 claims to take multiple infinitesimal generators $B={v_a}$ as input, creating a mismatch between the stated theoretical assumptions and the algorithmic design.
Moreover, Proposition 3.3 and Algorithm 1 define the dependent variable as one-dimensional, whereas Proposition 3.4 and the Reaction–Diffusion experiment involve multi-dimensional dependent variables.

- Minor issues.
(1) The paper should briefly explain the meaning of ordinary invariants when they are first introduced, to assist readers who are not familiar with the term.
(2) The subtitle Water Wave is inconsistent with the corresponding paragraph in the appendix and with the related figures and tables, which may cause confusion during reading.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 3
questions: - What precisely is new compared to DI-SINDy and EquivSINDy?

- Please add a baseline where a carefully curated physics dictionary for SINDy is compared against the invariant library.

- Please answer the weakness above.
soundness: 2
strengths: The authors present Propositions 3.3-3.4 that provide a way to integrate symmetry constraints to symbolic regression, and also a relaxed constraint mechanism, which is useful in practical applications.

The experiments consider realistic cases of noise, imperfect symmetries, and it is useful that their method integrates with different approaches because one can consider the symbolic regression method that is most suitable for the problem at hand, for example Weak-SINDy for noisy data. The success criteria are also clear. 

The paper is well-written and the method is well communicated.
summary: This paper presents an approach to enforcing symmetry constraints in symbolic PDE discovery. More specifically, the authors provide a Cramer’s rule type of constructive methodology to build higher order from lower order ones, a practical solution to real systems with imperfect symmetry, SI-relaxed, and ways to integrate the approach with multiple symbolic regression algorithms, more precisely sparse regression, genetic programming, and transformers, and demonstrates strong performance on three benchmarks.
weaknesses: The practical significance of the approach is limited by the following weaknesses:

- There exist a different paper on the same subject following a similar approach called “ Governing Equation Discovery from Data Based on Differential Invariants” which is neither discussed not compared against. The authors should provide an explicit technical comparison, e.g. explain how assumptions differ, how invariant extraction differs, and how the linear-subspace approach allows for the use of weak SINDy. Right now the methodological overlap might strike to some readers as significant. I would also recommend discussing the technical similarities and how this paper differs from EquivSINDy. This way the novelty can be made more clear. 

- The paper argues that considering symmetries shrinks the search space and improves accuracy and efficiency. For many PDEs, a domain-expert can define a small, physics-consistent library (knowing the experimental setup makes very accurate educated guesses possible) without considering any symmetry. The study should include a realistically curated library that is symmetry agnostic but contains a physically plausible set of terms to measure the benefit of considering invariant instead of physical process knowledge. In practice practitioners consider exactly these types of libraries. 

- Table 1 reports the complexity, the success probability, and the error between the predicted and the baseline PDEs. In several tasks, the reduction of the complexity is small, only a few variables, yet the narrative is strongly mentions complexity reduction as the main motivation. The authors need to carefully quantify how much the complexity reduction relies on invariants compared to other factors such as better choice of dictionary or better LHS choice. Also, define C in the main text and not in the Appendix because it is important. 

- Some invariant contain terms that are operations between derivatives. The authors augment baselines with PySINDY* libraries, but you need to be careful because if invariants can be expressed with degree-2 monomials in $\eta$ while in the jet variable you require a degree-3 then the comparison favors the invariant basis by construction. This is precisely why an “expert” defined dictionary needs to be considered, because one can consider many irrelevant terms to blow up the dictionary of SINDy that are not related to the physical process for which we try to discover a PDE.

- For the Boussinesq equation, a filter $| u_x | < 0.1$ is considered. That is a reasonable choice, but it introduces a data-selection bias that may increase or decrease the probability of finding the correct PDE. The authors should perform a study of how different values of the threshold, e.g. (0.01, 0.05, 0.1, 0.2, 0.3), affect the discovery and report the fraction of discarded samples. 


- The method assumes that the symmetry group is know and correct. This is acknowledged but the practical problem that discovering that full symmetry group from data is hard, and misspecification can result to the method braking down. The authors propose a relaxed approach for approximation, but there should a detailed study here. How wrong can a group selection be before the method collapses? Maybe they can run a sweep of different scenarios to provide some practical guidance. I also believe that the authors should consider symmetry braking due to forcing or BCs which is a major practical case. 

- The paper correctly states that any G-invariant PDE can be written purely in G-invariants. However this guarantees that the representation exists, it doesn’t guarantee that the representation is compact. In practice though an invariant algebra might require higher order generators or introduce expressions that a physics library wouldn’t so library size in the end might be comparable to the jet term library (see Table 1). The manuscript should discuss when the invariant set would be provably smaller and quantify the regime that reductions are substantial. 

- All three problems have nice symmetries, for example scaling and translations or SO(2). In this setting, some chosen invariants might be very close to the target structure, which makes the identification problem significantly easier than the jet and than the invariant in more complex regimens. 

- For me the invariant description is not interpretable meaning that if I know that a system contains diffusion or surface diffusion, I know which terms I need to include in the library, and how to prune it. Symmetry invariants are not as interpretable, unless you first start from the physics mechanisms and then define the invariants based on that. 

Overall I believe that this paper is mathematically elegant, but I am not convinced that it is useful in its current form.