Discovering Symbolic Differential Equations with Symmetry Invariants

We thank the reviewer for the constructive feedback. We address the question below:

A major advantage of the proposed approach is its compatibility with various existing algorithms, such as symbolic regression (SR) and SINDy. However, all the methods discussed in the paper rely on RMSE (with regularization) as the objective function. It would also be worthwhile to discuss and clarify in the paper whether the proposed approach can be extended to different objective functions, such as the variational objective used in D-CIPHER [1].

First, we would like to note that we have implemented our method on Weak SINDy, which uses a variational objective. Weak SINDy is used in all experiments in Sec. 4.4.

Nonetheless, this is a great suggestion, and we have conducted an additional experiment to implement our method on D-CIPHER and tested both the original D-CIPHER algorithm and our symmetry-enforced method on the Darcy flow dataset. The results are available in the table below.

We further explain the three methods in the table. The first method, D-CIPHER, uses the original implementation from the D-CIPHER official codebase. While their formulation of the extended differential operator in the paper is $a(\mathbf x) \partial_\alpha [h(\mathbf x, \mathbf u)]$, their implementation only addresses the operators in the form $\partial_\alpha [h(\mathbf x, \mathbf u)]$. This does not cover the terms like $xu_x$ and $yu_y$ in the Darcy flow equation. Therefore, in the second method of the table, D-CIPHER$^*$, we extended their code implementation by allowing function coefficients $a(\mathbf x)$ in the operator. Only in this way is the method theoretically capable of finding the correct functional form of Darcy flow. The final row represents our method, where we apply the procedure in Appx. C.3 (Line 778-809) to provide D-CIPHER with a set of rotationally invariant differential operators and fundamental rotational invariants, so that the spatial rotation symmetry is enforced in any discovered equation. More specifically, the inner optimization of D-CIPHER solves a 5-dimensional coefficient vector for 5 invariant operators in a constrained least-squares problem, and the outer optimization uses genetic programming to find a free-form expression of the 2 fundamental rotational invariants, $x^2 + y^2$ and $u$.

It can be seen from the table that our method can find the correct functional form, while D-CIPHER with the original variables and derivative operators cannot. The benefit of symmetry is even greater here for D-CIPHER than for other SR methods like SINDy, because D-CIPHER requires the user to specify both the function coefficient $a(\mathbf x)$ and the function to be differentiated $h(\mathbf x, \mathbf u)$ for an extended derivative (i.e., step 1 in their Figure 1). Such choices of functions can be largely arbitrary if no prior knowledge is available. On the other hand, our symmetry-based approach automatically selects this dictionary of differential functions.

In the revision, we will include this experiment and clarify that our approach can be extended to different objective functions, such as the variational objective.

We thank the reviewer for the constructive feedback. We address the questions/concerns below:

Presentational issues

We thank the reviewer for pointing out the presentational issues. In Sec. 2.1, we will elaborate on the intuition of prolonged group actions on jet spaces: to analyze the symmetry of PDEs, we must know how it transforms not just the variables, but also their derivatives. We will also include a visualization of the prolonged group action on a 2D space in the revision. We will also change the terminology in Appx. B, where “trivial” refers to “identical” transformation/function in that context.

As suggested by the reviewer, we will also reduce the paper’s reliance on cross-references to the appendix in the revision. We will make sure to include important implementation details and experiment results, such as those in Appx. C.1, C.2 and D.2, in the main text. We will remove the explicit link to equations in the Appendix in L222, and instead write out the full example.

L71: Elaborate on generalization to multiple variables and equations

Here, we mean our method can also be applied to discovering PDEs or systems of multiple PDEs containing more than one dependent variable ($q>1$). Several definitions/statements need to be extended in this case. For example, the infinitesimal generator in Eq. (3) becomes $\mathbf v = \xi^j \frac{\partial}{\partial x^j} + \phi^k \frac{\partial}{\partial u^k}$ to account for the multiple dimensions in $\mathbf u$. The computation of differential invariants still follows the same procedure; however, in Prop. 3.3, because the total space becomes $X \times U \simeq \mathbb R^p \times \mathbb R^q$, the number of lower-order invariants needed also becomes larger: $q_n = q\cdot \binom{p+n-1}{n}$. The rest of our method remains the same.

In Sec. 2, we restricted the setting to a single dependent variable to avoid introducing too much notational and conceptual complexity at once. This is used as the default setup in the paper. When the discussion of multiple dependent variables is needed, we explicitly mention such generalization, such as in Prop. 3.4. Again, this is a choice of how we present the paper – we could also make the most general statement at the beginning, possibly at the cost of increased presentational complexity and less clarity. We’d appreciate the reviewer’s thoughts on this.

Line 108/109: "we estimate the partial derivative terms" - how expensive is that? Would it be problematic for high-order differential operators?

In our experiments, we use the central finite difference method to estimate the partial derivatives (e.g., using numpy.gradient). The computational complexity is linear w.r.t. the dataset size and the differential order, so it is not expensive even for high-order differential operators. The problem for high-order derivative estimation is not computational cost, but robustness to noise. A small noise in the solution function would lead to a highly inaccurate estimation of high-order derivatives with the finite difference method. In that case, we use the variational objective proposed by WSINDy [1] to bypass the need for derivative estimation.

References

[1] Messenger, D. A. and Bortz, D. M. Weak sindy for partial differential equations. Journal of Computational Physics, 443:110525, 2021a.

Line 276: Why is $\eta(1,0)$ excluded?

Because $\eta(1,0)=\frac{u_x}{u_x}=1$ is a constant. Constant terms are already handled by symbolic regression algorithms. Thus, we do not add a constant function to the SINDy library or a constant feature to genetic programming.

Thank you for the constructive feedback. We address the questions below:

1. This method considers the scenario that symmetries are known while equations themselves are unknown. Is this a reasonable assumption?

Yes, this is a reasonable assumption. As we mentioned in Sec. 5, this assumption aligns with common practice – physicists often begin by hypothesizing the symmetries of a system and seek the most economical equations respecting those symmetries. Classic examples include deriving Navier-Stokes from rotational+translational symmetry and conservation laws, writing effective field-theory Lagrangians by enumerating all operators allowed by a symmetry group, etc. The assumption of symmetry is also commonly used in related works on equation discovery [1,2].

References

[1] Gurevich, Daniel R., et al. "Learning fluid physics from highly turbulent data using sparse physics-informed discovery of empirical relations (SPIDER)." Journal of Fluid Mechanics 996 (2024): A25.

[2] Messenger, Daniel A., Joshua W. Burby, and David M. Bortz. "Coarse-graining Hamiltonian systems using WSINDy." Scientific Reports 14.1 (2024): 14457.

2. In a real-world scenario, is there a situation where the equations are known, symmetries are unknown?

Such cases indeed exist. Some complex PDEs are written down before their full Lie-symmetry structure is classified. A separate line of work tackles this task of automated symmetry detection. We have cited those works in Sec. 5 and mentioned the integration of their methods and ours as future work.

While interesting, that scenario is orthogonal to our contribution: we show how given symmetry knowledge can be injected into any symbolic regression method to reduce the hypothesis space and guarantee physically valid outputs.

Note: We did not quite understand the transition from the previous question to this one. If the reviewer meant the reverse, please refer to the answer to Q1.

3. Isn’t this method just a reparameterization of the original variable by symmetry invariants?

No, our contributions are much more than that. Simply reparameterizing the variables is not enough because it is not compatible with some symbolic regression (SR) algorithms, such as GP which fits an explicit equation, and Weak SINDy which uses a variational objective. Thus, in addition to this general idea of using invariants:

• We developed Alg. 1 and Alg. 2 (in Appx. C.1) to implement this idea in different SR methods, with the additional procedure of LHS selection. As mentioned in Sec. 3.3, these adaptations are necessary because base SR methods have their own hypothesis spaces for equations, and we aim to implement the intersected hypothesis spaces (Fig. 2).
• For sparse regression (SINDy), we developed another technique to convert the symmetry condition (Thm. 3.2) to linear constraints on SINDy parameters. This allows us to impose symmetry on SINDy and its variants (e.g., Weak SINDy) easily. In this case, we do not reparameterize the variables, but use a derived theoretical result (Prop. 3.4) for easier implementation when it is applicable.
• In Sec. 3.4, we also described how to relax the symmetry constraint in our method. This extends the applicability of our method to systems with imperfect symmetry. The experiments showed that our “SI-relaxed” variant still outperforms an unconstrained baseline, showing that symmetry information, even if slightly misspecified, provides visible benefits.

4. Is this method applicable when the symmetries are more complicated, so that the symmetry invariants are difficult to compute? The paper mentioned a generic algorithm for computing symmetry invariants, but didn’t use it.

This is a misunderstanding. We did use the generic algorithm – see Appx. B.4, for example, for a full derivation for the rotation symmetry. And yes, the method is applicable in generic cases. The first step solves a first-order linear PDE (the characteristic equations for the Lie generators). This is algorithmic and implemented in standard computer-algebra systems such as sympy (pdsolve, for relatively simple systems) and Mathematica (DSolve). Then, the second step uses Eq. (6) from Prop. 3.3 to recursively construct higher-order functions. This reduces to symbolic differentiation and multiplication, again provided by off-the-shelf CAS libraries.

Regarding more complicated symmetry: although some classical PDEs, e.g. the 1D heat equation, admit large (even infinite-dimensional) Lie-symmetry algebras once the equation is known, our setting is the reverse: the equation is unknown and the symmetry is supplied a priori as an inductive bias. Hence, we typically posit only the symmetries that are evident from empirical observations or fundamental principles, such as spatial homogeneity (translations), isotropy (rotations), simple scalings, etc, not the full algebra that might emerge later. In other words, the groups we consider are often low-dimensional products of these elementary transformations, all readily handled by the method in Sec. 3.2. Moreover, should an unusually large group be encountered, the pipeline still functions; the only effect is a larger (but still finite) invariant basis.

5. Is this method applicable to more general, real-world equations?

Yes. To demonstrate this, we conduct another experiment on a simulated dataset for the Navier-Stokes equation. We follow the setup in [3], aiming to discover the vorticity equation $\omega_t=-\partial_x(\omega u)-\partial_y(\omega v)+0.01\Delta\omega$. The following table shows the results. Our method uses invariants of the rotation symmetry and has fewer variables and functional terms. As a result, it can discover a parsimonious equation matching the functional form of the ground truth.

We note that the result is based on our implementation of SINDy and our own dataset. The code and data for exactly reproducing the results in [3] are unavailable. Nonetheless, this result shows that our method using symmetry invariants is better than using original variables under the same backbone method implementation and data for equation discovery.

[3] Messenger, Daniel A., and David M. Bortz. "Weak SINDy for partial differential equations." Journal of Computational Physics 443 (2021): 110525.

6. The selected equations contain explicit invariant terms. Aren’t those well-chosen equations that work especially well with this method?

They are not deliberately chosen, but rather, those invariant terms arise naturally once the assumed symmetry holds. For example, the rotational invariance of Darcy flow implies isotropic permeability, and $x^2+y^2$ is the only ordinary function satisfying that property. These terms were not injected to “fit” our algorithm; they are the physically correct building blocks that any symmetry-respecting model would employ.

Also, our experiments tested beyond the perfect symmetry setting. In Sec. 4.4, we showed that even if there is symmetry breaking in the system and non-invariant terms in the equation, our method can still discover the correct equations via constraint relaxation. This further demonstrates that our method not only works on tailored problems.

7. The paper mentioned a generic algorithm for computing invariants, but in fact, it used precomputed invariants.

This is a misunderstanding. Please see our answer to Q5 for details. In short, the “precomputation” is exactly implementing the generic algorithm in Sec. 3.2.

8. Is it feasible to compute differential invariants using Prop. 3.3? E.g., there are determinant operations, which might cause numerical instability.

Yes, it is feasible. In Prop. 3.3, we compute the invariant functions instead of their evaluations on data points, which is purely symbolic and does not involve numerical issues.

We do note that evaluation of the precomputed invariant functions can cause numerical problems such as division by 0. However, handled properly, these numerical problems won’t affect the feasibility of the method. For example, we might encounter a rational function $f(x,u^{(n)}) / g(x, u^{(n)})$ as one of the invariants. In this case, we apply a filter $|g| > \epsilon$ for some small threshold $\epsilon$ to the dataset to discard large values that destabilize the regression. We will make sure we describe such procedures in the experimental details section in the revision.

9. How beneficial is the ability to handle non-sindy-format (not linearized) equations? Even when there are non-linearizable terms, series expansion may handle them.

The ability to handle “non-linearized” equations is critical. The series expansion of nonlinear terms introduces infinite power series and explodes the search space. While truncation is possible, we may still need higher-degree polynomials in order for the equation to match the data accurately. As a result, the regression is less efficient due to a larger function library, less accurate due to truncation errors, and the results are less interpretable due to the mix of expanded terms from nonlinear functions. Apart from our method, other related works, such as D-CIPHER [4] mentioned by other reviewers, also used genetic programming to handle non-SINDy-format equations and demonstrated its advantage over SINDy-based approaches.

[4] Kacprzyk, Krzysztof, Zhaozhi Qian, and Mihaela van der Schaar. "D-CIPHER: discovery of closed-form partial differential equations." Advances in Neural Information Processing Systems 36 (2023): 27609-27644.

We are happy to discuss further if there are remaining questions. If your concerns have been addressed, we kindly request that you adjust the evaluation accordingly.

I thank the authors for their responses, and some of my concerns (2, 4, 6, 7, 8, 9) are addressed. However, I have some concerns remaining about this paper.

1. Lie point symmetries are highly PDE-specific and can be quite nontrivial. Some are elementary—such as translations or rotations—whereas others are far more intricate. For example, the KdV equation $u_t + 6uu_x + u_{xxx} = 0$ admits the dynamical-scaling generator $3t\partial_t + x\partial_x - 2u\partial_u$ [1]. The cubic nonlinear Schrödinger equation $i\phi_t + \phi_{xx} + 2|\phi|^2\phi = 0$ has a “pseudo-conformal inversion” symmetry generated by
   $t^2 + tx \partial_x + \frac{i}{2}x^2 \phi \partial_{\phi} + t\phi \partial_{\phi}$ [2].

Because of this variety, identifying a Lie point symmetry is not as straightforward as, e.g., recognizing that the isometry group of $\mathbb{R}^{3}$ is $SO(3)$. Could you give concrete guidance—or name an explicit class of physical equations—for which your method applies? In other words, how large is the shaded region in your Venn diagram (Figure 2)?

2. The novelty of the proposed work is still uncertain.

(a) (SR) The algorithm 1 is a brute-force application (for LHS search) of SR algorithm using symmetry invariants.

(b) (Sindy) If I understood it correctly, the algorithm explained in the “sparse regression” section could be seen as parametrization of the symmetry invariants by linear combinations of library functions, by imposing linear constraints on W. This method seems novel, but its effectiveness seems to depend strongly on the particular library selected. Could you clarify how broadly applicable the method is and how sensitive it is to the choice of basis functions?

[1] Cantwell, Brian J. Introduction to Symmetry Analysis. Cambridge University Press, 2002.

[2] Devi, Preeti, and K. Singh. “Lie Symmetry Analysis of the Nonlinear Schrödinger Equation with Time Dependent Variable Coefficients.” International Journal of Applied and Computational Mathematics 7 (2021): 23.

Thank you for the follow-up. While we believe our previous response already discusses the two points you raised, we would like to elaborate a little further on those.

The method's focus on simple Lie point symmetries limits its applicability to more complex PDEs.

First, it is necessary to distinguish between the complexity of Lie point symmetries and the complexity of PDEs. Many complex PDEs also exhibit simple symmetries from first principles, e.g., the Navier-Stokes equation with spatial rotation symmetry. Thus, our method can be applied to complex PDEs in real world; see, for example, our initial answer to Q5 for the results on discovering the Navier-Stokes equation in vorticity form.

Then, if the reviewer is still concerned that our method mainly aims at simple symmetries, we would reiterate the explanations in our last response: these simple symmetries are ubiquitous, and their presence can be easily tested numerically, making them accessible inductive biases. Meanwhile, more complex symmetries are often PDE-specific, enumerating the possible forms of these symmetries is nontrivial, and it is difficult to obtain such knowledge before discovering the equation itself.

In terms of methodology, given any symmetry group, our pipeline of computing and using the invariant functions as regression features still works -- see our initial answer to Q4. But we intentionally choose not to consider more complex symmetries because such knowledge is not often available. In other words, this is not an intrinsic limitation of our method, but a limitation of available prior knowledge, which could possibly be addressed by incorporating symmetry discovery methods.

The novelty of the approach is still unclear, as it heavily depends on specific assumptions about the choice of basis functions. The method’s general applicability is not well-established.

We respectfully disagree. Consider our base method of using invariant functions instead of regular variables. It does not require additional assumptions about the choice of (SINDy) basis functions. It is general and can be applied to different SR algorithms and PDE datasets.

We understand that the reviewer has been questioning the novelty of this base method, and may not have referred to it when making the above statement. We would like to again highlight that this method of using invariants, alone, is a valuable and novel contribution recognized by all other reviewers.

Then, we understand that the above statement may specifically refer to the variant of our method that imposes linear constraints on SINDy parameters. We agree that it cannot be applied to any given symmetry group and SINDy basis functions. However, we have clearly established the conditions for it to work, which can be easily verified with a few symbolic calculations. We also show that common setups in related works on SINDy satisfy these conditions. And when these conditions are not met, we can always fall back to the base method, directly using invariant functions as features and targets for SINDy regression.

We hope the reviewer can take these clarifications into account when forming the final judgement. Thank you again for engaging in these discussions.

We thank the reviewer for the constructive feedback. We address the questions/concerns below:

The method relies on explicitly known symmetry groups without end-to-end learnable capability to extract symmetries from data. Thus, the application of this approach is limited.

Discovering symmetries and discovering equations from data are two distinctive tasks. Extracting symmetries from data was done by a separate line of work [1,2,3]. It is not the task we are solving in this paper.

We argue that using known symmetry groups is a reasonable problem setup. As we mentioned in Sec. 5, this aligns with common practice – physicists often begin by hypothesizing the symmetries of a system and seek the most economical equations respecting those symmetries. Classic examples include deriving Navier-Stokes from rotational and translational symmetry and conservation laws, or writing effective field-theory Lagrangians by enumerating all operators allowed by a symmetry group. The assumption of symmetry is also very common in related works on equation discovery [4,5].

In addition, the application of our approach is not limited to systems matching the given symmetry perfectly. We have also developed a method to relax the symmetry constraint in Sec. 3.4. This extends the applicability of our method to systems with imperfect symmetry. The experiments also showed that our “SI-relaxed” variant still outperforms an unconstrained baseline, demonstrating that symmetry information, even if slightly misspecified, provides a visible benefit.

References

[1] Yang, Jianke, et al. "Latent space symmetry discovery." arXiv preprint arXiv:2310.00105 (2023).

[2] Ko, Gyeonghoon, Hyunsu Kim, and Juho Lee. "Learning infinitesimal generators of continuous symmetries from data." Advances in Neural Information Processing Systems 37 (2024): 85973-86003.

[3] Hu, Lexiang, Yikang Li, and Zhouchen Lin. "Symmetry discovery for different data types." Neural Networks (2025): 107481.

[4] Gurevich, Daniel R., et al. "Learning fluid physics from highly turbulent data using sparse physics-informed discovery of empirical relations (SPIDER)." Journal of Fluid Mechanics 996 (2024): A25.

[5] Messenger, Daniel A., Joshua W. Burby, and David M. Bortz. "Coarse-graining Hamiltonian systems using WSINDy." Scientific Reports 14.1 (2024): 14457.

Comparison with D-CIPHER (Kacprzyk et al. (2023))

We agree with the reviewer that D-CIPHER is a relevant baseline and should be included. We ran D-CIPHER and also implemented our method upon D-CIPHER in an additional experiment. We tested the methods on the Darcy flow dataset in our paper and report the results in the table above.

First, we explain the three methods in the table. The first method, D-CIPHER, uses the original implementation from the D-CIPHER official codebase. While their formulation of the extended differential operator in the paper is $a(\mathbf x) \partial_\alpha [h(\mathbf x, \mathbf u)]$, their implementation only addresses the operators in the form $\partial_\alpha [h(\mathbf x, \mathbf u)]$. This does not cover the terms like $xu_x$ and $yu_y$ in the Darcy flow equation. Therefore, in the second method of the table, D-CIPHER$^*$, we extended their code implementation by allowing function coefficients $a(\mathbf x)$ in the operator. Only in this way is the method theoretically capable of finding the correct functional form of Darcy flow. The final row represents our method, where we apply the procedure in Appx. C.3 (Line 778-809) to provide D-CIPHER with a set of rotationally invariant differential operators and fundamental rotational invariants, so that the spatial rotation symmetry is enforced in any discovered equation. More specifically, the inner optimization of D-CIPHER solves a 5-dimensional coefficient vector for 5 invariant operators in a constrained least-squares problem, and the outer optimization uses genetic programming to find a free-form expression of the 2 fundamental rotational invariants, $x^2 + y^2$ and $u$.

It can be seen from the table that our method can find the correct functional form, while D-CIPHER with the original variables and derivative operators cannot. The benefit of symmetry is even greater here for D-CIPHER than for other SR methods like SINDy, because D-CIPHER requires the user to specify both the function coefficient $a(\mathbf x)$ and the function to be differentiated $h(\mathbf x, \mathbf u)$ for an extended derivative (i.e., step 1 in their Figure 1). Such choices of functions can be largely arbitrary if no prior knowledge is available. On the other hand, our symmetry-based approach automatically selects this dictionary of differential functions.

Even if symmetry is enforced by using symmetry invariants in the SR algorithms, the resulting symbolic equations can still be very complex and not necessarily physically meaningful, especially when dealing with noisy data and imperfect symmetry.

While it may be true that our method can produce complex equations in the presence of imperfect symmetry or just generally challenging setups, this type of complexity is the result of the underlying system, not a failure mode of our algorithm.

Our method aims to find the equation, ideally more parsimonious, that can explain the data well. But if the data itself comes from a complex system, then our method would unavoidably need to include more terms in the symbolic equation. E.g., for imperfect symmetry, we have developed a method to relax the symmetry constraint and discover symmetry-breaking components (Sec. 3.4) and showed how it worked experimentally in Sec. 4.4. The discovered equations indeed contain more terms, but they fit and explain the data better than simpler equations. Also, the additional terms are still physically meaningful, as they correspond to different ways symmetry can be broken in that reaction-diffusion system.

The reviewer also mentioned noisy data. We believe this is a somewhat different challenge – the underlying system can be simple, but the data can be noisy due to inaccurate measurements. However, in this scenario, our symmetry-based approach still works well with weak-form SR algorithms that improve robustness to noise, including Weak SINDy (Sec. 4.4) and D-CIPHER (as described above).

How can the proposed method avoid complicated solutions with excessively high orders that lose physical meaning?

The maximal derivative order is a hyperparameter of our proposed method. We only compute invariant functions up to this order and use them as input features to symbolic regression. We agree that excessively high-order derivatives lose physical meaning and make regression more difficult by introducing irrelevant features. We will include this discussion in the revision.

If the reviewer meant something different by “high orders”, we would appreciate clarification.

Please let us know if our response has addressed your concerns. We are happy to discuss further if there are remaining questions. If your concerns have been addressed, we kindly request that you adjust the evaluation accordingly.