Discovering Symbolic Partial Differential Equation by Abductive Learning

We thank Reviewer B51F for the insightful review, the careful examination of our paper, and the appreciation of the value and innovation of our work. Below, we answer your questions in detail, hoping to address your concerns.

---

Q1. Clarification on constants in the grammar and additive coefficients.

A1. Thanks for your question. In this work, the only constants are the linear combination coefficients in Equation (1). These constants are not generated when abducing candidate terms; instead, they are introduced during FEMTrain as trainable parameters. To introduce a constant source, it suffices to add an additional column of ones and assign it a trainable coefficient. We will revise Equation (1) to explicitly include the additive coefficient for clarity.

---

Q2. Concern regarding the re-write of complex terms.

A2. Thanks for your question. The growth rate in the number of terms caused by re-writing is indeed an interesting question. While deriving a formal bound is difficult due to the lack of a clear definition for a compact term, we can offer an intuitive explanation.

First, it is more appropriate to consider rewriting sets of complex expressions, since a PDE typically contains multiple terms. Consider a set of complex expressions A and a set of simpler, canonical expressions B that can subsume A. Since the full space of complex expressions contains exponentially more elements than the space of canonical terms, it is clear that the size of B does not necessarily grow exponentially with the size of A; otherwise the full space of simpler canonical terms would be larger. Moreover, we generally prefer an increase in the number of terms rather than an increase in the number of symbols per term. This is because the search space grows exponentially with the number of symbols, whereas increasing the number of terms does not expand the full search space of MSS, which is desirable for our approach. We will add this discussion to the next version of our paper.

---

Q3. Clarification on the FEM-based residual and operator equivalence in Theorem 1.

A3. Thank you for your insightful comment. We recognize that our use of the term 'weak formulation' was misleading, for which we sincerely apologize. You are correct that our FEM-based residual does not contain the boundary terms expected from a traditional weak formulation.

To clarify, our FEM-based residual is more accurately described as a Galerkin projection of the strong-form residual. We compute derivative terms (e.g., $u_{xx}$) directly using automatic differentiation on a neural network, and then project the resulting strong-form PDE residual onto the FEM basis functions. Because this procedure does not involve integration by parts, boundary integrals do not arise in this part of our loss function. The boundary values of the solution are enforced via the MSE component of the consistency measure defined in Equation (2). We will revise the manuscript to reflect this precise description of our approach.

Regarding your other concern on Theorem 1, you are correct that its application as PDE operators relies on certain assumptions. In our framework, the crucial assumption is that the solution is sufficiently smooth for the required derivatives to exist. This smoothness is provided by the neural network, which acts as a continuous and differentiable representation of the solution field. The derivatives are then evaluated via automatic differentiation. We will clarify in our revision that the practical application of Theorem 1 rests on this smoothness being effectively captured by the neural network backbone.

Thank you again for your valuable feedback, which will significantly improve the clarity and accuracy of our paper. Your insight is also very enlightening for us. We consider the development of a solver capable of automatically performing symbolic weak form computations, as well as methods for the adaptive selection of finite element basis functions, to be important directions for our future work.

---

Q4. Clarification on the computational cost and numerical stability of the LevelMargin step.

A4. Thanks for your question. The projection step in MSS is a standard least-squares problem, and the dimension of the subspace projected onto each time does not exceed the initial sparsity. Since there is a closed-form solution, the computation is very efficient. In our implementation, we use PyTorch’s method for computing matrix pseudo-inverse, which robustly handles colinear terms. We did not observe any numerical instability in our experiments. As for the condition number of the mass matrix M, we do not need to be concerned about M being ill-conditioned. This is because we are free to choose the mesh for the FEM-based residual in the consistency measure. "Certainly, the effect of the condition number of M on discovery performance is a valuable question. We consider this an important point for our subsequent research.

---

Q5. Clarification on the proof of Lemma 5.

A5. Thank you very much for the detailed review! You clearly pointed out that when $v$ is not orthogonal to $Y$, the left side of the inequality will change. However, please note that this does not cause the inequality to fail. We omitted some details in the original proof, which we now provide in full.

Proof:

Let $v$ be a vector that is not orthogonal to the space $\text{span}(Y)$. Decompose $v$ into two orthogonal components: $v = v\_Y + v\_{Y^\perp}$, where $v\_Y = \text{Proj}_{\text{span}(Y)}(v)$ is the component of $v$ in $\text{span}(Y)$, and $v\_{Y^\perp}$ is the component orthogonal to $\text{span}(Y)$.

Since $\text{span}(Y \cup \{v\}) = \text{span}(Y \cup \{v\_{Y^\perp}\})$, the RHS remains unchanged:

$$\\|\text{Proj}\_{\text{span}(Y \cup \{v\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(Y)}(u)\\|\_2 = \\|\text{Proj}\_{\text{span}(Y \cup \{v\_{Y^\perp}\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(Y)}(u)\\|\_2.$$

For the LHS, consider the relationship between $v_Y$ and $\text{span}(X)$:

1. Case 1: $v_Y \in \text{span}(X)$
   If $v_Y$ lies in $\text{span}(X)$, then $\text{span}(X \cup \{v\}) = \text{span}(X \cup \{v_{Y^\perp}\})$. The inequality becomes: $$\\|\text{Proj}\_{\text{span}(X \cup \{v\_{Y^\perp}\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(X)}(u)\\|\_2 \ge \\|\text{Proj}\_{\text{span}(Y \cup \{v\_{Y^\perp}\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(Y)}(u)\\|\_2.$$
2. Case 2: $v_Y \notin \text{span}(X)$
   Decompose $v_Y$ further as $v_Y = v\_{Y,X} + v\_{Y,X^\perp}$, where $v\_{Y,X} = \text{Proj}_{\text{span}(X)}(v_Y)$ and $v\_{Y,X^\perp}$ is orthogonal to $\text{span}(X)$. Then: $$v = v\_{Y^\perp} + v\_{Y,X} + v\_{Y,X^\perp}.$$ The subspace $\text{span}(X \cup \{v\})$ simplifies to $\text{span}(X \cup \{v\_{Y^\perp}, v\_{Y,X^\perp}\})$. According to the Pythagorean theorem, projecting onto a larger subspace increases the norm and therefore the following inequality holds: $$\\|\text{Proj}\_{\text{span}(X \cup \{v\})}(u)\\|\_2 \ge \\|\text{Proj}\_{\text{span}(X \cup \{v\_{Y^\perp}\})}(u)\\|\_2.$$

From Lemma 5:

$$\\|\text{Proj}\_{\text{span}(X \cup \{v\_{Y^\perp}\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(X)}(u)\\|\_2 \ge \\|\text{Proj}\_{\text{span}(Y \cup \{v\_{Y^\perp}\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(Y)}(u)\\|\_2.$$

Combining results, we have:

$$\\|\text{Proj}\_{\text{span}(X \cup \{v\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(X)}(u)\\|\_2 \ge \\|\text{Proj}\_{\text{span}(Y \cup \{v\})}(u)\\|\_2 - \\|\text{Proj}\_{\text{span}(Y)}(u)\\|\_2.$$

This confirms the inequality in all cases.

---

Q6. Clarification on the comparison with D-CIPHER.

A6. Thanks for your question. D-CIPHER is a SOTA PDE discovery method. One of its characteristics is to integrate measurements of the same PDE under multiple initial and boundary conditions, which likely contributes to its robustness to noise. However, since our benchmarks do not involve multiple datasets for one PDE, we set the number of datasets to 1 for D-CIPHER to align the experimental setup and repeated each experiment 5 times, reporting the best result. To address your concerns about the experimental results, we have re-run D-CIPHER using the original paper’s hyperparameter settings (10 datasets, 10 trials) and report the results using their metric, RMSE, as shown in the table below.

The experimental results align well with those reported in Figure 3 of the original paper. We will include these updated results and provide a clearer description of the experimental setup in the next version of the paper.

---

Q7. Clarification on the candidate term selection.

A7. We use the same knowledge base to abduce candidate terms for all tasks, with the only difference being that we state different dependent and independent variables for each task.

---

Q8. Clarification on the discovery result under 5% noise.

A8. Thank you for your detailed review of our paper! This non-monotonic result likely arises from our consistency optimization (CO) stage for two reasons:

1. Hyperparameter $\lambda$: We used the same $\lambda$ across different noise levels. It's plausible that the chosen $\lambda$ is slightly better suited for a 5% noise measurement.
2. Noise as Regularization: A small amount of noise can act as a regularizer, preventing overfitting and improving the model's generalization. Better generalization is critical for discovering a generalizable PDE, sometimes more so than perfectly fitting clean data.

Notably, the physical field fitting error shown in our paper's Figure 3 remains monotonic with noise, which supports the overall validity of our experiments.

Thank you for redesigning the margin. This removes the orthogonality issue nicely.

Now you just need to do a few minor things to finalize the revision i) restate Lemma 5/ Theorem 2 under the new definition and provide a short proof that making $F \subseteq H$ can only reduce $|| \text{Res} ||_2$, ii) provide a short clarification on how the fixed library $H$ is chosen and how sensitive is the method to the $0.05$ threshold, iii) comment on whether OMP-style guarantees apply to your method or not.

With these additions all my concerns about the soundness of the theoretical part will be fully resolved and I will increase score to acceptance.

We thank Reviewer dxqN for detailed review, helpful comments, and appreciation on the practicality and innovation of our work. Below, we answer your questions in detail.

---

Q1. "In Table 2, the authors compare the performance of ABL-PDE and its variants with the baselines PDERidge and DL-PDE++. I have a few questions and comments: ... Do you ignore the coefficients of redundant terms and only use the correct terms when computing the error?"

A1. Thanks for your helpful question. We will revise our paper to clarify these points.

1. Exactly. We used the same pre-trained neural network and the same regression method to ensure a fair comparison.
2. The standard deviation presented here is with respect to our method's hyperparameter $\lambda$, which is not applicable to other methods. As mentioned in Section 6.2, results for DL-PDE++ were obtained through careful tuning of its hyperparameters. When we are unable to tune its hyperparameters to discover the correct equation, we ensure it includes all ground truth terms and minimize redundant ones. For our method, we try to include relevant terms as much as possible at the outset by setting a very small margin threshold ($1 \times 10^{-3}$), and redundant terms will then be eliminated later through consistency optimization.
3. Thanks for your comments. We list the number of redundant terms below and will add this to the next version of our paper (empty cells indicate no redundant terms):

4. Exactly. When computing the error, we ignore the coefficients of redundant terms and only use the correct terms.

---

Q2. "Could the authors compare their method and its variants with more recent and representative baselines? ... It would be helpful to include comparisons with more modern methods such as Symbolic Physics Learner (Sun et al., ICLR 2023) or SGA-PDE (Chen et al., Physical Review Research 2022)."

A2. Thanks for your suggestion. These are indeed two highly influential works. We extensively discussed SGA-PDE in our paper, treating it as a representative of asynchronous methods. While Symbolic Physics Learner (SPL) doesn't directly involve PDEs, its core idea of using MCTS to handle large search spaces is orthogonal to our approach and offers significant inspiration for our future work. We will add a discussion of SPL in the next version of our paper. Since the code for SPL only covers ODEs rather than PDEs, applying it to our benchmark is non-trivial. Therefore, we invested effort to adapt ABL-PDE to ODEs and conducted experiments on their benchmark, Lorenz system. Similarly, because the code for SGA-PDE does not support the discovery of PDE systems, while 2 out of 3 of our benchmarks are systems of PDEs, we also applied ABL-PDE to the benchmarks used by SGA-PDE. We summarize the results in the table below. The results of SGA-PDE and SPL are taken from their respective papers. The results of our methods are obtained with the same initial sparsity $10$ and margin threshold $1 \times 10^{-3}$ as in the paper. Due to the tight rebuttal schedule, we did not select $\lambda$ based on the fitting error of the physical field on the validation set. Instead, we simply used a trick to rescale the supervised loss and the physics loss to similar magnitudes for all the experiments: $\text{loss} = \text{loss}_1 + \frac{\text{loss}_2}{\text{loss}_2 / \text{loss}_1 + 1e-8}$.

References:

[1] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). Physical Review Research, 2022.

[2] Symbolic Physics Learner: Discovering governing equations via Monte Carlo tree search. ICLR, 2023.

---

Q3. "Since the overall optimization process alternates between MSS and FEMTrain, can the authors provide a proof or argument that this process always converges?"

A3. Thanks for your question. A formal convergence proof is indeed challenging, since the overall process alternates between two complex subproblems: structure identification (MSS) and parameter/coefficient optimization (FEMTrain). Proving convergence for such alternating schemes is notoriously difficult, especially when the subproblems are not convex, which is the case here. For future work, we will consider analyzing the convergence of our method under reasonable assumptions.

---

Q4 & Q5. "The caption of Figure 3 mentions that 5 different values of lambda ... Does this mean five separate configurations with $\lambda$ set to 0.5, 0.75, 1.0, 1.25, and 1.5 respectively? The wording is a bit unclear."

A4 & A5. Thanks for your question. We will clarify this point in the next version of our paper. The following table shows $\lambda$ for each experiment:

We thank Reviewer htDg for insightful review and appreciation on the originality of our work. Below, we answer your questions in detail.

---

Q1. "can you add add a synthetic example where the underlying PDE has unknown nonlinear parameter functions to see how the current framework fails or could be extended."

A1. Thanks for your question. As mentioned in the conclusion, the limitation of ABL-PDE lies in its inability to handle unknown parameters within nonlinear terms. For example, $\theta$ in $\frac{1}{\theta \cdot u + v}$ is a typical case. It should be noted that introducing such terms into the hypothesis space of our knowledge base is straightforward; the main challenge lies in computing the consistency measure, particularly in the Monotone Subset Selection (MSS) step. In our work, the minimal computational unit of MSS involves solving for the linear combination coefficients $\xi$ given a set of expression terms. This is a least squares problem with a closed-form solution, allowing for efficient computation. However, if unknown parameters are introduced within nonlinear terms, the minimal computational unit of MSS becomes an optimization problem that may lack a closed-form solution and can even be non-convex, resulting in significantly increased computational cost. This is precisely why we mentioned in our paper the necessity to incorporate advanced techniques for joint parameter optimization for future work.

---

Q2. "Though standard PDEs were used, the inclusion of a more realistic engineering problem (e.g., turbulent flow, reaction-diffusion with unknown kinetics) could strengthen the practical impact."

A2. Thanks for your suggestion. We highly agree that testing our method in more complex scenarios would better demonstrate its practical applicability. However, fitting unknown and complex kinetics often involves learning parameters within nonlinear terms. As we mentioned in A1, learning these parameters significantly increases computational cost and requires more advanced optimization techniques, which we consider as future work.

---

Q3. "If the initial pre-training fails (very sparse or chaotic data), the abductive loop may not recover."

A3. Thanks for your question. We acknowledge that such scenarios present significant challenges, and this is indeed an inherent difficulty in the broader problem of scientific discovery from data.

In situations with very sparse data, where a dense grid is unavailable for traditional numerical differentiation, employing neural networks for function approximation is a widely adopted technique [1-2]. Our bidirectional loop acts as an enhancement to these data-driven approximation methods. It allows us to iteratively refine the model's understanding of the physical field by leveraging discovered physical insights, even if initial approximations are coarse due to data scarcity.

For chaotic systems, methods with stronger inductive biases or prior knowledge—such as Gaussian Processes, Fourier series, or Chebyshev polynomials—might indeed offer better performance than a generic neural network. Fortunately, our framework is highly modular and flexible. As long as these alternative models incorporate learnable parameters (e.g., length scales and signal variance for Gaussian Processes), they can seamlessly substitute the "machine learning" part depicted in Figure 2 of our paper. This adaptability allows our overall system to leverage the strengths of different modeling paradigms, tailoring the data approximation component to the specific characteristics of the physical system.

References:

[1] DLGA-PDE: Discovery of PDEs with incomplete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 2020.

[2] Deep-learning based discovery of partial differential equations in integral form from sparse and noisy data. Journal of Computational Physics, 2021.

We thank Reviewer DouY for detailed review and suggestions on the clarity of our paper. We will revise the paper accordingly. Below, we answer your questions in detail.

---

Q1. "Why is it necessary to have a consistency maximization process?"

A1. Thanks for your question and suggestion. The intuition behind the consistency maximization process is to bridge the gap between the continuous, differentiable machine learning component and the discrete, non-differentiable logical reasoning component, as illustrated in Figure 2 of the main text.

We can train the neural network with gradient descent to approximate the physical field, however, we cannot calculate the gradient of the equation structure with respect to the neural network output. Thus, we requires a calculable criterion to construct a connection between the equation structure and the neural network, such that the change of neural network parameters can quantitavely affect the equation structure. In simple words, the consistency measure makes the equation structure an implicit function of the neural network parameters, and we optimize the neural network and determine the equation structure jointly by maximizing the consistency measure.

---

Q2. "In the contribution, the authors mention bidirectional loop, however symbolic regressions method like genetic programming also has this kind of iterative loop structure, how is the authors’ method different from that?"

A2. Thanks for your question. We highly agree with you that many symbolic regression methods, such as those based on genetic programming, also employ iterative loop structures. However, the fundamental distinction of our method lies in the nature and purpose of its bidirectional loop.

Specifically, our bidirectional loop operates between neural network training and equation structure identification. This contrasts sharply with the iterative loop structures found in other methods, which are confined solely to the realm of equation structure identification. As we mentioned in lines 51-53 of our paper, other approaches ignore the potential of using discovered physical information to enhance the model's ability to fit the physical field.

This distinction applies to genetic programming-based methods like DLGA-PDE [1], DL-PDE [2], and SGA-PDE [3], as well as sparse regression-based methods like WeakSINDy [4] and D-CIPHER [5]. To the best of our knowledge, our method is the first to establish a mutually beneficial loop between model training and equation structure identification.

References:

[1] DLGA-PDE: Discovery of PDEs with incomplete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 2020.

[2] Deep-learning based discovery of partial differential equations in integral form from sparse and noisy data. Journal of Computational Physics, 2021.

[3] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). Physical Review Research, 2022.

[4] Weak SINDy for partial differential equations. Journal of Computational Physics, 2021.

[5] D-CIPHER: discovery of closed-form partial differential equations. NeurIPS, 2023.

---

Q3. "What is the input/output of the neural network used for derivative estimation? There are many methods for approximating derivatives and the authors mention their limitations in the introduction. How does having a neural network mitigate those issues?"

A3. The input of the neural network is the spatial-temporal coordinate and the output is the predicted value of the physical field at that coordinate.

We mitigate issues of derivative estimation through our bidirectional enhancement loop (i.e., consistency optimization). The core idea here is to leverage the discovered physical information to guide and improve the neural network's training. Even if this physical information is initially imprecise, the resulting improvements in the network's performance can, in turn, help refine the physical information itself, creating a mutually beneficial and iterative enhancement.

Experimental results demonstrate the effectiveness of consistency optimization in improving the accuracy of derivative estimation. As discussed in our analysis of the Burgers' Equation experiment, the shock wave in the solution can only be captured by an extremely dense mesh that is impractical in most real-world applications. This limitation leads to suboptimal derivative estimation performance for the pre-trained network. Even when the form of the equation is known (as in PDERidge, Table 2), the coefficient estimation still have significant errors. In contrast, our method achieves a 5x to 10x reduction in L2 relative error through consistency optimization.

---

Q4 & Q5. "Can the authors give an example of how abduction works as indicated in section 4.1, step 2. It is currently hard to understand how abduction relates the number of operators to expand the library of candidates. ... Table 1 integrity constraints are hard to follow at present, despite the authors giving details. Can you please elaborate more on one constraint, maybe C1 as to why this will give a unique order."

A4 & A5. Thanks for your question. Due to space limitations, we use a simplified example—multiplication expressions with two dependent variables (u, v)—to detail the processes of abduction and the role of integrity constraints.

First, we present a vanilla Prolog program without operator number limit and integrity constraints.

% Logical facts
unit_expr(u).
unit_expr(v).

% Logical rules
multiplication_expr(X) :- unit_expr(X).
multiplication_expr(mul(X, Y)) :- unit_expr(X), multiplication_expr(Y).

The semantic of the program is intuitive: a multiplication_expr is either a unit_expr or a mul with two arguments-the first being a unit_expr, and the second being a multiplication_expr. When querying ?- multiplication_expr(X)., Prolog attempts to satisfy this goal by finding a value for X that fits the definition of a multiplication_expr. It finds answers using Selective Linear Definite clause resolution (SLD resolution), in other words, through recursion and backtracking. We explain the process in detail below.

1. First Rule: Base Cases

   • Prolog first tries this rule. It looks for simple "unit" expressions.
   • It finds unit_expr(u) and reports X = u.
   • Then it finds unit_expr(v) and reports X = v. These are the simplest, non-nested expressions.

2. Second Rule: Recursive Cases

   • After the simple cases, Prolog moves to the second rule to build more complex expressions.
   • This rule says a multiplication_expr can be mul(something, something_else) where the first part is a unit_expr (like u or v), and the second part is another multiplication_expr.
   • How it works:
     • It picks a unit_expr for the first part.
     • Then, for the second part, it recursively calls multiplication_expr again. This means it goes back to step 1 and 2 to find another expression for the second part.
     • This recursion allows it to build nested structures like mul(u, u), then mul(u, v), mul(u, mul(u, u)), and so on.
     • Backtracking: If Prolog hits a dead end or you ask for more solutions (by typing ;), it backtracks. This means it goes back to the last choice it made and tries a different option to find more combinations and deeper nested structures.

Operator Number Limit

One issue with the above program is that it will generate endlessly nested structures like mul(u, mul(u, mul(u, mul(u, ...)))), and another subtle issue is that the number of v will never be greater than 1. By introducing operator number limit, we solve both issues.

unit_expr(u).
unit_expr(v).

multiplication_expr(X, Op_limit) :- 0 =< Op_limit, unit_expr(X).
multiplication_expr(mul(X, Y), Op_limit) :- 1 =< Op_limit, unit_expr(X), 
                                            multiplication_expr(Y, Op_limit - 1).

The only difference is that we add a new argument Op_limit to the multiplication_expr rule, which is the number of operators allowed in the expression. When querying ?- multiplication_expr(X, 1)., only u, v, mul(u, u), and mul(u, v), mul(v, u), mul(v, v) will be generated.

Integrity Constraints

Since we do not want to both generate mul(u, v) and mul(v, u), we introduce the following integrity constraints to the second rule.

multiplication_expr(mul(u, v), Op_limit) :- 1 =< Op_limit, unit_expr(u), 
                                            multiplication_expr(v, Op_limit - 1),
                                            \+ integrity_constraint(mul(X, Y)).
integrity_constraint(mul(X, Y)) :- Y = mul(Z, _) -> X @> Z ; X @> Y.

The semantics of \+ is "not," meaning that for a valid mul(X, Y), the boolean value of integrity_constraint(mul(X, Y)) should be false.
The definition of integrity_constraint(mul(X, Y)) uses an if-then-else structure (Condition -> Then ; Else), which evaluates to true if the former term is greater than the latter term according to the meaning of @>, a built-in predicate in Prolog.
By introducing this integrity constraint, we can prevent the generation of terms like mul(v, u).

For a more comprehensive example of our logic-based knowledge base implementation, please refer to Appendix B.