Discovering Symbolic Partial Differential Equation by Abductive Learning
Authors: En-Hao Gao
Cunjing Ge
Yuan Jiang
Zhi-Hua Zhou
Venue: NeurIPS 2025 poster ()
Keywords: Abductive Learning
AI4Science
Primary area: machine_learning_for_sciences
TL;DR: We propose a dual-driven (data and knowledge) approach for symbolic PDE discovery that exhibits superior performance in terms of noise robustness and hyperparameter stability.

Abstract: Discovering symbolic Partial Differential Equation (PDE) from data is one of the most promising directions of modern scientific discovery.
    Effectively constructing an expressive yet concise hypothesis space and accurately evaluating expression values, however, remain challenging due to the exponential explosion with the spatial dimension and the noise in the measurements.
    To address these challenges, we propose the ABL-PDE approach that employs the Abductive Learning (ABL) framework to discover symbolic PDEs.
    By introducing a First-Order Logic (FOL) knowledge base, ABL-PDE can represent various PDEs, significantly constraining the hypothesis space without sacrificing expressive power, while also facilitating the incorporation of problem-specific knowledge.
    The proposed consistency optimization process establishes a synergistic interaction between the knowledge base and the neural network learning module, achieving robust structure identification, accurate coefficient estimation, and enhanced stability against hyperparameter variation.
    Experimental results on three benchmarks across different noise levels demonstrate the effectiveness of our approach in PDE discovery.

Discovering Symbolic Partial Differential Equation by Abductive Learning
En-Hao Gao, Cunjing Ge, Yuan Jiang, Zhi-Hua Zhou
National Key Laboratory for Novel Software Technology. Nanjing University, China School of Artificial Intelligence, Nanjing University, China gaoeh@lamda anjmma.n...... edu. cn, {gecunjing, jiangyuan, zhouzh}anju. edu.
Abstract
Discovering symbolic Partial Differential Equation (PDE) from data is one of the most promising directions modern scientific discovery. Effectively constructing an expressive yet concise hypothesis space and accurately evaluating expression values, however, remain challenging due the exponential explosion with the spatial dimension and the noise in the measurements. To address these challenges, we propose the ABL-PDE approach that employs the Abductive Learning (ABL) framework to discover symbolic PDEs. By introducing a First-Order Logic (FOL) knowledge base, ABL-PDE can represent various PDEs, significantly constraiming the hypothesis space without sacrificing expressive power, while also facilitating the incorporation of problem-specifie knowledge. The proposed consistency opti- mization process establishes a synergistic interaction between the knowledge base and the neural network learning module, achieving robust structure identification, accurate coefficient estimation, and enhanced stability against hyperparameter variation. Experimental results on three benchmarks across different noise levels demonstrate the effectiveness of our approach in PDE discovery.
Introduction
Partial Differential Equations (PDEs) serve as fundamental mathematical tools for describing vast array of physical phenomena in science and engineering. They effectively capture the intricate space and/or evolve over time. Complemented
relationships between how quantities change by both analytical techniques and modern numerical solvers, PDEs provide interpretable models that support reliable prediction and control. Due to their versatility, PDEs widely employed across a range of applications, including airfoil design in aerodynamics meteorology and pricing analysis in quantitative finance
14
weather prediction in
Recently, physics-inftee informedd machine learning has demonstrated remarkable success in approximating the behavior of complex physical systems [3] 22], which significantly changes the way we represent and leverage PDEs. Note that single PDE can be adapted to various scenarios by simply modifying its boundary conditions and domain geometries; however, PDEs in explicit form remain essential in applications where high reliability is paramount. Therefore, discovering symbolic governing equations from data becomes of the most attractive research directions in scientific discovery [4.
18 have been made to PDE discovery task. Most of these methods consist of two primary components: derivative-integra calculator and candidate term library, as illustrated in Figure The calculator evaluates symbolic expressions based on either direct derivative estimation or weak formulation-base integral computation. The library constructs the hypothesis space for the target equation, with each term representing function of the variables of interest. The underlying PDE is typically assumed to be linear combination of these terms and is identified through symbolic regression techniques.
Significant efforts
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

A grand challenge of PDE discovery lies in constructing a comprehensive candidate term library. On one hand, this library must be sufficiently expressive to capture pat- terns within data; On the other hand, it should remain compact to enable an effi- cient PDE structure identification. Exist- ing studies attempt to achieve the trade-off by constraining the expressiveness using a concise library [T] [6], or emphasizing
8,u
31
22
T
8ee = 0,u
Measurements
Discovered PDE
y 0,30
Derivative Integral
Candidate
Calculator
Library
Figure 1: An illustration of PDE discovery framework.
the expressiveness at the expense of using
a redundant library g 26 27). Accurately evaluating the expression values is also challenging. Because
only access discrete measurements, derivatives integrals have be estimated using approximate methods, the accuracy of which can be significantly affected there is no sufficient data or there is noise. Many studies attempt address this issue using techniques such as ensemble learning [12] and weak formulation [II], while this paper attempts to exploit the discovered physical information help enhance the calculator's is capability, which has rarely been exploited before.
In this paper, we propose the ABL-PDE (Abductive Learning for PDE discovery) approach to tackle the aforementioned challenges under the ABductive Learning (ABL) framework [28] (2)] ABL is a powerful paradigm that integrates data-driven machine learning with knowledge-driver logical reasoning in balanced mutual promotion loop while maintaining the expressive power of both. Building on this framework, we design two components for ABL.-DEF a neural network-hased learning module for derivative estimation and First-order Logic (FOL.)-based reasoning module for candidate term generation, along with a consistency optimization process to bridge the two components. Our approach balances the expressiveness and compactness of the hypothesis space in PDE discovery and enables reciprocal enhancement between the calculator and the library.
Our main contributions are summarized as follows:
We propose a novel PDE discovery approach based on the ABL framework, transforming the previous unidirectional discovery pipeline into bidirectional enhancement loop.
We introduce a FOL knowledge base capable of representing a wide range of PDEs, sub- stantially reducing redundancy and facilitating the incorporation of task aask-specific knowledge.
We design consistency optimization process to bridge the two components, which sig- nificantly enhances the robustness of structure identification, the accuracy of coefficient estimation, and the stability against hyperparameter variations in PDE discovery.
We conduct extensive experiments on various PDE discovery tasks across different noise levels. The
experimental results demonstrate the effectiveness of our approach.
2
Related Work
According to whether the equation structure and its coefficients are determined synchronously or not, PDE discovery methods can be roughly categorized into two classes.
The first class comprises synchronous methods, most of which rely on sparse regression techniques. A pioneering approach is SINDy [S which was extended by PDE-FIND 23 from ODEs to PDEs. Weak SINDY [11 further enhances noise robustness using weak formulation. Apart from explicit sparse regression, several synchronous methods such as PDE-Net [18 [9] and Bayesian hidden physics models M optimize derivatives and coefficients end-to-end. Although generally exhibiting superior performance in fitting physical ficlds, they struggle to distinguish noise from small-coefficient terms. D-CIPHER [66 extends the boundary of PDE discovery to learn parameters within nonlinear terms; this comes with a challenging optimization process, making it struggle to learn coupled PDE systems. Generally speaking, synchronous methods rely on predefined, typically polynomial libraries that ensure term uniqueness, though such compact library may lack expressiveness when encountering complex structures. Actually, they may even find it challenging to identify fractions.
The second class comprises asynchronous methods, often based on evolutionary algorithms, that do not need be constrained by predefined libraries g 12 L5 25). These methods first determine the equation's form by discrete optimization and then estimate the equation's coefficients. Although

---

such approaches incur higher computational costs due to discrete optimization, they mitigate the detrimental impact interference terms on coefficient estimation while facilitating the incorporation
of more flexible equation representation methods, such as binary trees [
15
25
Note there are
great progresses in establishing theoretical foundation for evolutionary algorithms [0], such that
many of these methods are not purely heuristics any more. These methods, however, often generate an excessive number redundant and invalid PDEs within the hypothesis space. By contrast, through introducing FOL knowledge base, our method achieves better balance between expressiveness and compactness in PDE representation, while also facilitating the incorporation of prior knowledge.
3 Preliminaries
In this section, we present the problem setting and introduce the abductive learning framework.
3.1 Problem Setting
In this paper, we consider general dynamical systems governed by PDEs of the following form:
d.u= { (1) R)), x)) (u, 2)), ..... .).....
p (w, 2), E:
(1)
where u [0, x Ω + R is the solution function with time horizon T and bounded spatial domain c. R4, t and æ denote the temporal and spatial coordinate, each f is function of u and æ, and their linear combination coefficient which is always sparse in practice. The data available is set of triplets {t, Z4, ü(t,, | representing m temporal-spatic coordinates, with measurements that may contain noise. Each ( x) collectively forms the candidate term library that will be gradually expanded during the execution of our approach. We call f(+) (u, x) an expression term.
Definition (Expression Term). An expression term is recursively formed by the following rule:
(expr) :== (idpv) (dpv) (uspp)) (expr) (expr) (op) (expr)
(op) ::= +
//
Here, (-) denotes a nonterminal symbol, ::= indicates definition, and signifies alternation. The nonterminal symbols (dpv) and (idpv) represent dependent and independent variables, respectively, Combining Equation (T) with expression terms, suffices to represent a wide range of fundamental PDEs, such the Diffusion equation, the Schrödinger equation, the Navier-Stokes equation, etc.
3.2 Abductive Learning
Abductive Learning (ABL) 28 29 is powerful paradigm that integrates data-driven induction from machine learning with knowledge-driven deduction from logical reasoning. A key challenge this integration is that the discrete deductive process is non-differentiable ABL overcomes this by using abductive reasoning to propagate the discrepancy between the logical deductions and data labels back to the machine learning model, which then serves the supervision signal for its training.
Abductive reasoning. also known as abduction, is fundamental mode of logical inference that seeks the best explanation for observations. A common approach to implement this process is through Abductive Logic Programming (ALP) [7], where an abductive logic program is formally defined by triplet (P, A, IC'. P is the background knowledge encoded as logic program, A is set of abducible predicates representing hypotheses, and /C' a set of first-order formulae called integrity constraints. Abductive explanation for observation O is set A of ground atoms on abducible predicates A, such that the following conditions are met (where -= denotes logical entailment):
(1) PUAF O,IC
(2) PUA A is consistent
4 Our Approach
In this section, we first present an overview of our approach, ABL-PDE (ABductive Learning for PDE discovery), and then detail its components and the proposed consistency optimization process.
3

---

Data
Candidate Terms
Operator Number
Orr u
2. Abduction
2
Train
3
I.. FEMTrain
4.MSS
NN
Oz
KBB
Evaluate
82 arx
i
O
Machine Learning
Logical Reasoning)
Expression Value
Figure 2: An illustration of ABL-PDE's framework.
4.1 Overview of ABL-PDE
As illustrated in Figure 2 ABL-PDE comprises two modules:
Machine Learning. This module employs a neural network to fit the physical field and calculates derivatives through automatic differentiation.
Logical Reasoning. This module incorporates an abductive logic program as the knowledge base and is responsible for generating candidate expression terms.
The discovery process commences with pretraining the neural network, during which no physical information is available (Step 1). Following this initial stage, an iterative deepening search is performed by progressively increasing the number of operators in the candidate terms. In each iteration, ABL-PDE first expands the candidate term library through abduction, utilizing the current number of operators observations (Step 2). After this expansion, the abduced expression terms transformed into numerical values (Step 3). Next, the proposed Monotone Subset Selection (MSS) algorithm employed to identify terms with nonzero coefficients (Step 4). Subsequently, ABL-PDE estimates their coefficients and jointly optimizes the neural network through Finite Element Method- based Training (FEMTrain) (Step 5). Steps and together constitute the consistency optimization process, which is repeated until the discovered equation stabilizes and the coefficients converge.
4.2 Logic-based Knowledge Base
This part details our knowledge base KB, an abductive logic program (P, P,A,,1 A, IC).
Background Knowledge P. A straightforward way is to adopt the syntax rule in Definition aas the background knowledge P. Since P U A must be consistent, P effectively defines the boundary of the search space by constraining A to be valid expression term. This space. however, grows exponentially with the number of operators, posing a significant challenge. We observe that this vast space contains substantial redundancy, many terms are mathematically equivalent despite differing in their syntactic representation. To mitigate this issue, we introduce canonical form term, which prune the search space by eliminating these equivalent representations.
Definition 2 (Canonical Form Term). An expression term is called in canonical form if and only if it can be recursively generated by the following rules:
(cfterm) ::= (div) (mul) (add) (sub) (unit)
(div)
(mul) (add) | (mul) >
(sub)
(mul) : (unit) (mul) (unit)
(add) (mul) + (add) (mul)
(sub)
(add) (add)
(unit)
(dpv) (idpv) (par)
(par)
O(app)) (dpv) i (ssapy) (par)

---

Table 1: General-purpor (CI, C2, C3) and domain-specific (C4, C5, C6) integrity constraints.
Integrity Constraint
Description
CI add(X, L) L aYd Y 2) A X Y.
Unique addition order
C2 + mul (X, L) AL- = mul() (Y, 2) A X
Unique multiplication order
C3 f- sub(X,Y) intersect(X, Y).
Disjoint subtraction operands
C4 getPower(mul(X, Y); u, N) N n.
Power limit
C5 < getOrder(par X, M) A M > m.
Derivative order limit
C6 getDependent (aar =Z=
Continuity equation constraint
The notation has the same meaning as in Definition These rules impose structure constraints on expression terms derivable from each nonterminal symbol. We illustrate the constraint on partial derivatives ((ar)) in Example 1 for detailed explanations of other rules, please refer to Appendix
Example 1. Given two dependent variables and The expression term 0 (u x is not in canonical form does not satisfy the rule of (par) nor any other rules. However, 0 (uxv) = - which is the linear combination of canonical form terms 0r and u x 0,v. In Section prove that every expression term can be represented by a linear combination canonical form terms.
In our implementation, P logic program directly translated from the rules Definition Integrated with the operator number tracking capability in the deductive reasoning process, P can leverage abductive reasoning to generate canonical form terms.
Abducible Predicates A. In this work, the design of abducible predicates is not focus. We simply introduce an auxiliary predicate to only allow abducing terms derivable from (mul) and (div) Since we consider linear combinations of expression terms, this does not compromise the expressiveness.
Integrity Constraints IC. Integrity constraints enable to declaratively state various constraints on the structure of expression terms. Table [ presents several examples of such constraints: the upper section lists general-purpose constraints, while the lower section includes domain-specific ones. The logical symbols are used as follows: 'r"' denotes implication (with the left-hand side False omitted), 'A' represents conjunction, and ">" indicates symbolic order, which built-in predicate in most logic programming languages. CI and C2 in Table [ state that an addition or multiplication sequence invalid if the former of two adjacent subterms symbolically greater than the latter. C3 ensures that the two operands of subtraction do not intersect. C4 and C5 respectively constrain the power and the derivative order with respect to x term. C6 is pertinent the continuity equation. Specifically, in our experiment on the Navier-Stoker equation, knowing that 2D incompressible flow satisfies the continuity equation (-ō,,u Öy 1) allows us eliminate terms involving partial derivatives of with respect to y, since these terms can be systematically replaced by partial derivatives of with respect (e.g., By - -Ug, Byy = -ury). Please note that only the general- purposse constraints are included in KB to ensure a fair comparison with other methods.
4.3 Consistency Optimization
In this part, we first introduce the consistency measure and then present MSS and FEMTrain.
Consistency Measure. Given the operator number limit L and the dataset D, the consistency measure Con (NNo, KB) of neural network NNo and knowledge base KB is defined as follows:
- ( min Res &, P) MSE(NN,, D))
s.t. FC U AxR,O,
0=1
where 9 denotes the neural network parameters, & is the coefficient vector as in Eq. ( F a set of terms obtained from the abductive reasoning KB for each operator number 0 < L, and is hyperparameter utilized to balance Res and MSE.
Pleasse refer to Appendix B for pedagogical example illustrating the program and the abduction process.
5

---

Let F = (f x) h) the physics informed phosics-informed residual is defined as:
NF & f(() (t
Res( N,NNo x))
Let D = {(t,, Ri ait;, x.))) as in Section
3.1 the mean square error is defined as:
m
MSE(NNg, D) m
(NNo(t;, x1) ū(t,, xx))
The minimal residual in Eq. decreases monotonically as the size of F increases. Thus, directly maximizing the consistency measure will lead to a trivial solution where UAa.o Since the number of nonzero terms cannot be determined in advance, selecting
upper bound of [F] is insufficient to balance sparsity and consistency. Thus we introduce margin and k-level margin. Definition 3 (Margin). Let H be set of expression terms, and F C H. The margin of F with respect to H is defined as:
Margin (F) |Res* (NNo, F)1l, |Ress (NNe, H)ll
(5)
|Res (NN,, H)ll2
where Res" (NNo. is the minimal residual in Eq. with respect &.
Definition 4 (k-Level Margin). Let H be set of expression terms, k positive integer, and ||| denote set cardinality. The k-level margin of H is defined as
LevelMargin
min F'CH, |
Margin,;("),
(6)
Intuitively, the margin a subset F' quantifies how self-sufficient it in approximating the dynamics 0,NNg. A small margin suggests that the terms within F are sufficient, while the terms outside it are negligible. Accordingly. the k-level margin measures the sufficiency of the best possible k-term model, indicating whether parsimonious, k-term equation shown in Section the level margin decreases monotonically as & increases, which enables us to adaptively control the expression term number by applying margin threshold.
accurately represent the system. As
Monotone Subset Selection. Under a given margin threshold, evaluating the con- sistency measure Con (NNg, KB) with re- spect to the neural network parameters 0 time-consuming, as the complexity the residual minimization subproblem in Eq. (2) is no less than that of the sparse regression under cardinality constraint, which is generally NP-hard [20]. To reduce computational costs, propose Mono- tone Subset Selection (MSS), a two-stage algorithm featuring an initial stage of rapid, large-scale coarse selection succeeded by a precise, small-scale refined search.
Algorithm MSS
Input: Neural network NN Term library termLib, Ini- tial sparsity ini initSp, Margin threshold magThd Output: Minimal residual set under BagThd-mmrghh minSet
I: Initialization: minSet + 0
libVal < Eval NN, termLib)
3: subLib <- POSS(libVal, termLib, initSp)
left +- 1, right < | |subLib
5: while 1eft right
6:
mid + []eft toveft+rrii right)/2]
7:
tmpMinSet, levMag + Enum mid) if levMag < magThd then
8:
9:
minSet tmpMinSet, right mid
Algorithm] I presents the details of MSS. 10:
else
the first stage (cf. Line B). MSS evaluates 11: left mid
expression terms in the library and calls
Pareto Optimization for Subset Selection (POSS) 21 30 find reduced term library subLib whose size should be no greater than the initial sparsity. POSS a non-deterministic subset selection algorithm designed to minimize the residual in Eq. (2) under cardinality constraint. In expectation, is guaranteed to find a reasonably good solution within a polynomial number of queries to the objective function. This property makes it well-suited for the first stage, which focuses on rapidly identifying relevant terms rather than performing costly residual minimization. More details about
SSinee the underlying PDE is expected to hold at every point within the domain, we are free to select the set of points for calculating the residual, which has been omitted for brevity.

---

POSS are provided in Appendix [ In the second stage (cf. Line MSS employs binary search find the smallest subset size for which the level margin below the given threshold. Each time, enumerates the subset of size mid to identify the one with the minimal residual and calculates the level margin, which is tractable within the reduced library size.
Finite Element Method-basee Training. Once the set of active expression terms F' in Eq. (2) is identified via MSS for a given neural network parameters 0, the next step of consistency optimization is to jointly optimize the network parameters and coefficient vector & to maximize the consistency measure. A standard approach for this sub-problem to adopt the Physics-Informed Informed Neural Networks (PINN) paradigm [22]. However, the efficacy of the PINN approach is often hindered by its sensitivity to hyperparameters. such as the loss weighting term À in Eq. (2) [3]. To enhance the stability of model training and coefficient estimation, we move beyond pointwise collocation and instead enforce the physics-informer residual
weighted-average sense across the domain. This achieved [24] of the residual from Eq. (3) onto the space spanned by Finite
through Galerkin projection
Element Method (FEM) basis functions, leading to the following objective function:
(O Σ (NN (t, x))
where M is the mass matrix constructed from the inner product of FEM basis funetions. Please refer to Appendix for detailed definitions and derivations. The overall consistency optimization process alternates between MSS and FEMTrain until the F stabilizes and the coefficients & converge.
5 Theoretical Analysis
In this section, we first demonstrate that our logic-based knowledge base, KB, retains full expressive power, and then establish the soundness of MSS. Proofs are available Appendix E
We denote the abduction result of KB for given operator number as H4, Thus, the set of all possible abduction results is given by H U H4. Besides, we refer to the set of all expression terms as S.
Definition 5 (Subsume). A finite set of expression terms F is said to subsume expression term there exists an expression term s' such that:
s' can be derived from through finite sequence of algebraic operations and differentiation. s' can be symbolically expressed as linear combination of expression terms in F, i... EEEF of f, where a are scalars.
As shown in Example 1 { u u,W d.u) subsumes On X v).
Theorem 1. For every E S, there exists a finite subset H' C H such that H' subsumes
According to the PDE formulation in Eq. (T. it suffices to consider the expressiveness of a candidate term library up to a linear combination. Theorem states that every expression term in has mathematically equivalent form that is linear combination of elements from H, which directly implies that the expressive power KB is equivalent to that of the full expression term space S.
The soundness of MSS depends on the criterion utilized its binary search, whose validity is ensured by the following monotonicity theorem.
Theorem 2. For expression term set F and integers 1 S kn lo S [F), it holds that:
(kx))Margin,,k:)) 2 LevelMargin,
6 Experiments
In this section, we present empirical studies to answer the following questions:
1. How does ABL-PDE perform against contenders in PDE discovery?
2. How robust is ABL-PDE against hyperparameter selection?
3. How does each component of ABL-PDE contribute its performance?
7

---

Table 2: PDE discovery comparison results. The evaluation metric is the sum of L2RE (%) on all coefficients. Number with a dagger (f) indicates the discovered equation has redundant terms. Results of ABL-PDE (w/o FEM) and ABL-PDE are the mean and standard deviation of 5 experiments with different À that used balance the two sub-metrics in Eq. D
We bold the best results in each task.
Hurgers' Eeguation
Schrüdinger Equuntion
Naiver- Stokes Equation
Method
5%
10%
0
5%
10%
5%
10%
PDERidge
22.68
23.86
24.54
10.34
10.26
17.13
16.58
40.92
DLE+
22.68
23.86
24.54
10.34
10.26
9.99
57.251
37.96'
98.52!
ABI-PDE (wo (O)
22.68
23.86
24.54
42.16"
44.677
48.76'
17.99
40.92
ABL-PD (wlo PKM)
13,09 + 2.81
14.65 + 2.67
15.04 + 3.87
3.05 + 0.94
3.27 + 1.09
2.87 + 1.07
30,88 + 5.24
26.47
5.00
34.28 + 7.50
ABL-DE
2,18
1,06
2.79 0.13
4.12 + 0.34
1,03 + 0.10
1,12 0,06
1,08
0.29
18.02 + 0.08
15.22
0,08
21.27 + 1.72
6.1 Experimental Setup
We conduct experiments on three types of discovery tasks: (i) oie-dimensional PDE: Burgers' Equation, (ii) one-dimensional PDE system: Schrödinger Equation, (ii) two-dimensional PDE system: Navier-Stokes Equation. Tasks (ii) and (iii) each involve two coupled equations. We use data from public datasets (2). incorporating 5% and 10% Gaussian noise to simulate inaccurate measurements as (27] We use the L2 Relative Error (I.2RF, %) as in [1] to measure the quality of coefficient estimation and physical field approximation. More details are specified in Appendix F Code is available at: https: / github..bbb... httpshhh..ooomeeceeeeeooe.eeeeee /wwucttiveeeerrrnnnn/aaaaa ABL-
6.2 Compared Methods
We compare ABL-PDE with four methods: 1) PDERidge: We build this strong baseline by assuming a known PDE structure and applying ridge regression to estimate coefficients. This represents the best possible performance for most methods that focus on enhancing structure identification [8] [0 [6].
2) DL-PDE++: DL-PDE (7 is a representative method among synchronous approaches that discussed in Section We enhance it by applying a normalization trick and tuning the sparsity threshold to retain all correct terms while minimizing redundancy. 3) ABL-PDE (w/o Con): utilizes the pre-trained neural network and MSS algorithm discover PDE, without the subsequent consistency optimization. 4) ABL-PDE (w/o FEM): It employs the residual in Eq. 3 instead of the FEM residual in Eq. 7 during the consistency optimization. To ensure a fair comparison, we use the same pre-trained neural network as in ABL-PDE to compute derivatives for the first two methods.
To verify the effectiveness of our basic knowledge base KB in reducing redundancy, we compare it with two baselines: 1) Full Space: The full space of expression terms as defined in Definition [ 2) Canonical Form: It contains all expression terms in the canonical form, which lacks abducible predicates and integrity constraints compared to KB.
6.3 PDE Discovery Comparison
In this part, we seek answer the first two questions. Table presents the results of structure identification and coefficient estimation. Number with a dagger (f) indicates the discovered equation has redundant terms (See Appendix for for the specific form). The performance of DL-PDE++ and ABL-PDE (w/o CO) is equivalent to PDERidge when the discovered structure matches the ground truth, due their shared neural network for derivative computation. As shown in the table, ABL-PDE correctly identified the PDE structure across all three tasks under varying noise levels, achieving the best coefficient estimation in 8 out of 9 cases and ranking second only to the structure-knows baseline, PDERidge, in the remaining one case. Notably, this performance was achieved utilizing the same initial sparsity of 10 and the same margin threshold of 0.05 for all three tasks. ABL-PDE has only one remaining hyperparameter, which will be discussed in Section 6.4 In our experiments, despite the simple structure of Burgers' equation, remains challenging for coefficient estimation. The inferior performance of compared methods is primarily attributed to the formation of sharp shock waves in the solution, which can only be captured by an extremely dense mesh that is impractical in most real-world applications. ABL-PDE effectively mitigates this issue via consistency optimization, which utilizes the discovered physical information, even if it initially incorrect, to improve the derivative estimation, leading to 5x to 10x reduction in L2RE.
8

---

ABL-PDE (w/o CD)
ABL- (w/o FEM)
ABL-
Test Enor (8)
Test Emor (%)
4.5
3.0
4.0
3.5
3.0
1.0
5
10
10
Noise Level (%)
Noise Level (%)
Noise Level (%)
(a) Burgers' Equation
(b) Schrödinger Equation
(c) Naiver-Stokes Equation
Figure 3: Test error comparison results. The evaluation metric is the sum of L2RE on each component of the predicted physical field. Results of ABL-PDE (w/ FEM) and ABL-PDE show the mean and standard deviation of experiments under different used in Eq. (2). The area the shadowed circle represents the magnitude of the standard deviation. Detailed results are provided in Appendix G.2
Beyond its robustness noise, ABL-PDE is also resilient hyperparameter selection. As shown in T the discovered equation of ABL-PDE (w/o CO) contains redundant terms in several cases, This redundancy stems from both the neural network's 's imperfect approximation to the physical field and, more critically, from a misalignment between the margin threshold and the data. Nevertheless, ABL-PDE successfully eliminates these terms through consistency optimization, demonstrating its notable robustness to sub-optimal hyperparameter choices. A by-product of our method is the improvement in the neural network's ' generalizability. Figure presents the L2RE of model prediction on the noise-free test data. Compared with the best performance achievable by supervised training alone (cf. ABL-PDE (w/o CO)), ABL-PDE consistently outperforms it across all tasks. Initially, the neural network performs suboptimally, and the candidate term library is too expansive. Consistency optimization progressively clarifies the ambiguity within the library, leading to a refined equation that, in turn, provides physical insights to improve neural network.
6.4 Ablation Study
Influence of FEMTrain. We investigate the effect of FEMTrain by comparing the performance of method with and without the utilization of FEM residual. Selecting hyperparameter, e.g., À in Eq. (2). to balance the physical information and the supervision signal has long been an existing issue with physics-infomne informed machine learning [13]. Therefore, each task with hyperparameters: [0.5A, 0.75X, X, 1.25A, 1.5A]. The mean and standard deviation of coefficient errors are presented in Table 2 and the results of test errors visualized in Figure 3 In addition to the superior performance of ABL-PDE in coefficient estimation and test error, its stability hyperparameter variations particularly noteworthy. In Figure we use the area of the shadowed circle to reflect the magnitude of the standard deviation. Compared to ABL-PDE (w/o FEM), the standard deviation of ABL-PDE is reduced by approximately one order of magnitude, This enhanced stability is crucial for practical scientific discovery, as finding an effective hyperparameter is much easier than finding the optimal one.
conduct five experiments for
Influence of KB. We use the one-dimensional PDE with showcase to demonstrate
N Term Number)
Full Space
single dependent variable the effectiveness of KB, i.e., the abductive logic program based on the canonical form, in constraining the hypothesis space. As the full space of expression terms is enormous, we use dynamic programming to compute its size, and the sizes of Canonical Form and ABL-PDE are counted by the implemented logic programs, respectively. Figure 4 illus- trates the logarithmic growth curve of the candidate term number as the operator number limit increases. Compared to the full space of expression terms, the proposed canon- ical form reduces the problem size to a computationally feasible scale. Nevertheless, KB still demonstrates a reduction by constant factor in the exponent, making it tractable for POSS. As terms with 10 operators can construct highly complex equations, our method can be applied to most scenarios.
Exceods range
Canonical Form
10
Operator Number Limit
Figure 4: Logarithmic growth curve of the candidate term number as the opera- tor number limit increases.

---

7 Conclusion
In this paper, we propose a novel PDE discovery method based on the abductive learning framework, aiming at leveraging the logic-based knowledge base and the neural network in a reciprocal way. Specifically, we introduce the ABL-PDE (ABductive Learning for PDE discovery) approach, which alternatively refines the knowledge base by utilizing neural network to access data and enhances the neural network's s performance by integrating the physical information provided by the knowledge base. Experimental results demonstrate the robustness of our approach in PDE structure identification, the accuracy in coefficient estimation, and the stability against hyperparameter variations. Furthermore, ABL-PDE is a general-purpont approach with sufficient flexibility in implementation, e.g., the neural network can adopt any advanced architecture, and the knowledge base can incorporate domain- specific information by declaring restriction to learning linear combination coefficients, which prevents from handling unknown parameters within nonlinear terms. In future work, we will try to construct parameterized candidate term library and incorporate advanced techniques for joint parameter optimization.
integrity constraints. The limitation of ABL-PDE lies in its
Acknowledgements
This research was supported by the Noncommunicable Chronic Diseases-National Science and Technology Major Project (2024ZD053I 1802, the National Natural Science Foundation of China (62202218), and the Jiangsu Science Foundation Leading-edge Technology Program (BK20232003).
References
[1] Steven Atkinson. Bayesian hidden physics models: uncertainty quantification for discovery of nonlinear partial differential operators from data. arXiv:2006.0422 04228. 2020.
[2] Reza Akbarian Bafghi and Maziar Raissi. PINNs-Torch: Enhancing speed and usability of physies-informed neural networks with pytorch. In The Symbiosis of Deep Learning and Differential Equations III, 2023.
[3] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate medium-range global weather forecasting with 3d neural networks. Nature, 61977970)5533-533, 2023.
[4] Steven L Brunton and Nathan Kutz. Promising directions of machine learning for partial differential equations. Nature Computational Science, 4(7):483-494 2024.
[5] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. National Academy of Sciences, 1 113(15):3932-399 2016.
[6] Martin Burger, Luis Caffarelli, and Peter A Markowich. Partial differential equation models in the socio-economic sciences. Royal Society A, 372(20099:211344066 2014.
[7] Alejandro Carderera, Sebastian Pokutta, Christof Schütte, and Martin Weiser. Cindy: Conditional gradient-based identification of non-linear dynamics-noise-roust recovery. 0rXiv:2101.02633 2021.
[8] Yuntian Chen, Yingtao Luo, Qiang Liu, Hao Xu, and Dongxiao Zhang. Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). Physical Review Research, 4(2)0223174, 2022.
[9] Jean Coiffier. Fundamentals of Numerical Weather Prediction. Cambridge University Press, 2011.
[10] Miles Cranmer. Interpretable machine learning for science with PySR and SymbolicRegres- sion.jl. arXiv:2305.01582 2023.
[11] A. Messenger Daniel and M. Bortz David. Weak SINDY for partial differential equations. Journal of Computational Physics, 443:110525,, 2021.
10

---

[12] Urban Fasel, J Nathan Kutz, Bingni W Brunton, and Steven L Brunton. Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control. Royal Society A, 447882260):202199044 2022.
[13] Zhongkai Hao, Songming Liu, Yichi Zhang, Chengyang Ying, Yao Feng, Hang Su, and Jun Zhu. Physics-informee machine learning: A survey on problems, methods and applications. arXiv:221100866, 2022.
[14] Raymond M Hicks and Preston A Henne. Wing design by numerical optimization. Journal of Aircraft, 15(7):407-412, 1978.
[15] Zhongyi Jiang, Chunmei Wang, and Haizhao Yang. Finite expression methods for discovering physical laws from data. arXiv preprint arXiv:2305.08342, 2023.
[16] Krzysztof Kacprzyk, Zhaozhi Qian, and Mihaela van der Schaar. D-CIPHER: discovery of closed-form partial differential equations. NeurlPS, 2023.
[17] A. C. Kakas, R. A. Kowalski. and F. Toni. Abductive logic programming. Journal of Logic and Computation, 2(6):719-770) 1992.
[18] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-Net: Learning pdes from data. In ICML, 2018.
[19] Zichao Long, Yiping Lu, and Bin Dong. PDE-Net 2.0: Learning pdes from data with numeric-symbolic hybrid deep network. Journal of Computational Physics, 399: 2019.
[20] Balas Kausik Natarajan. Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2):227-234, 1995.
[21] Chao Qian, Yang Yu, and Zhi-Hua Zhou. Subset selection by pareto optimization. In NeurIPS, 2015.
[22] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics- informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Compulational Physics, 378:686-707, 2019.
Physics-informed
[23] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances. ((4):e602261 2017.
[24] Gilbert Strang and George J. Fix. An Analysis of the Finite Element Method. Wellesley- Cambridge Press, 2nd edition, 2008
[25] Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discovering governing equations via monte carlo tree search. In ICLR, 2023.
[26] Hao Xu, Haibin Chang, and Dongxiao Zhang. DLGA-PDE: Discovery of PDEs with incom- plete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 418: 109584, 2020.
[27] Hao Xu, Dongxiao Zhang, and Nanzhe Wang. Deep-learning based discovery of partial differential equations in integral form from sparse and noisy data. Journal of Computational Physics, 445: 10592, 2021.
[28] Zhi-Hua Zhou. Abductive learning: towards bridging machine learning and logical reasoning. Science China Information Sciences, 62(7):76101, 2019.
[29] Zhi-Hua Zhou and Yu-Xuan Huang. Abductive learning. In Pascal Hitzler and Md. Kamruzza- man Sarker, editors, Neuro- -huro-Symboliccooo Artificial Intelligence: The State of the Art, pages 353-369. IOS Press, Amsterdam, 2022.
[30] Zhi-Hua Zhou, Yang Yu, and Chao Qian. Evolutionary Learning: Advances in Theories and Algorithms. Springer, 2019.
11

---

NeurIPS Paper Checklist
The checklist is designed encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No] or [NA]. .
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justificatin right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation.
While "[Yes] " is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] "provided proper justification is given (c.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used"). "[No] or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
general, answering
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS Paper Checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only the provided macros for your answers.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims the abstract and introduction accurately represent the paper's contributions and scope. They are supported by the thcoretical analysis and experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
12

---

Justification:
We discusses the limitations of the work in the conclusion section.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specifiitiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-texi system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: We provide the full set of assumptions and complete (and correct) proof for each theoretical result in the paper.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We provide detailed experimental settings Appendix and a detailed README file reproduce the results the supplementary material.
Guidelines:
13

---

The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We use public datasets and provide code with a detailed README file to reproduce the results in the supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //mips.cc/publcc//uut hww.b//iccguiddess/ddeeubbbbib
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
14

---

At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide detailed experimental settings in Appendix F and a detailed README file reproduce the results the supplementary material. We also analyze the influence of hyperparameters on the performance of the proposed method in the main text.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification:
We report error bars for the proposed method in the main text.
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We provide the detailed experimental settings in Appendix which contains the discussion of the compute resources used.
15

---

Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and confirm that the research conducted in the paper conforms
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper presents work whose goal is to advance the field of Machine Learning. There many potential societal consequences of work, none which we feel must be specifically highlighted here.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (c.g., deployment of technologics that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
16

---

Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We properly credit the creators of the assets and mention the license and terms of use explicitly.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version asset used and, possible, include URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, wwwpsswieeesseeeeeeeeee thcode com datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
We will release the code upon acceptance.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
17

---

At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
We do not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification:
We do not involve research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing. editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification:
The paper does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non standard components.
Please refer to our LLM policy http::///eerrrp.cc/ooofffffefee htt.:////ccc//oeeee/e///////2/ for what should or should not be described.
18

---

A Detailed Explanation of the Canonical Form
The canonical form of expression terms defined by the BNF grammar presented in Definition 2 This grammar uses standard notation: (-) for nonterminal symbols, ::= for production rules (indicating that the nonterminal on the left can be replaced by the expression on the right), and signifies alternation between possible forms. The notation used here is consistent with that in Definition [ In this grammar, nonterminal symbols define the valid structure for terms corresponding to specific operators or types: (cfterm) for any canonical form term, (div) for division (/), (mul) for multiplication (x), (add) for addition (+), (sub) for subtraction (-), (par) for partial differentiation (0), (unit) for base multiplicative/derrvatii
units, and (var), (dpv), (idpv) represent variables.
The production rules in Definition 2 define the set of expression terms that are in canonical form by imposing syntactic constraints on their structure. We explain the meaning of these rules as follows:
The top-level rule, (cfterm) :== (div) (mul) (add) (sub) (unit), specifies that any expression in canonical form must be derivable from one of these five main nonterminals or be a base unit.
The rule for division, (div) := (mul) / (add) (mul) / (sub), defines the valid structure for division terms (/). It requires the numerator be term derivable from (mul) and the denominator to be a term derivable from either (add) or (sub). For example, a term corresponding to (u + v)/(u t v) would be considered invalid under this rule because its numerator, (u + v). is derivable from (add), not (mul).
The recursive rule for multiplication, (mul) ::= (unit) x (mul) (unit), defines the structure of multiplication terms (x). The base case is any term derivable from (unit). The recursive step requires a term derivable from (unit) on the left of "x" and term derivable from (mul) on the right. This imposes specific structure on multiplication chains; for instance, a product of three terms like u x u x w can only be derived in the form (unit) (mul), where the second term derivable from (mul) is u w. This enforces the structure u x (V x w), preventing syntactically distinct but equivalent forms like (u x v) x w.
The recursive rule for addition, (add) :: (mul) + (add) (mul), similarly defines the structure of addition terms (+). It allows a base case of a term derivable from (mul), and a recursive step (mul) + (add). This ensures a specific structure for addition chains; for example, an expression like L + u x w (where X v is derivable from (mul) can only be derived as (mul) + (add). where the first term derivable from (mul) is and the term derivable from (add) is v + This enforces the structure t + (u v + w), preventing forms like (t + u x v) + w.
The rule for subtraction, (sub) := (add) - (add), defines subtraction terms (-), requiring both the left and right operands to be derivable from (add). This enforces a specific structure; for example, u + w (where u + v and w are derivable from (add)) can only be derived as (u + v) w, not (w v).
The rule (unit) := (dpv) (idpv) (par) defines terms derivable from (unit). These are referred as unit expressions and serve as the base elements for multiplication and differentiation arguments.
The rule for partial differentiation, (par) := (spp)) (dpv) 0: (par), defines partial derivative terms (ő). It requires the derivative operator be applied to a term derivable from either (dpv) or (par), and differentiation is always with respect to a term derivable from (idpv) (an independent variable). This imposes strong constraints, for instance, forbidding derivatives of products or sums as single terms, such as 0z (u x v). which would be invalid because u is neither derivable from (dpv) nor (par).
B A Pedagogical Example of Logic-based Knowledge Base
We provide a pedagogical example to elucidate the logic-based knowledge base and the abduction process discussed the main text. As a pedagogical example, this knowledge base considers three
basic operators (+, -,/). ,/), and the full implementation (reasoning/expp) of the supplementary material. By querying ?- write_ta to_ ttltttoffflll file. SWI-Prolop we can generate all expression terms with up to given number of operators (2 in this code). Abduced expressions
pl) is provided in the code
19

---

all in the form similar to mul (u, mul v)). A single integrity constraint is employed to ensure that multiplication expressions are a canonical order. illustrating how various human knowledge can be declaratively stated through such constraints. The declarative nature simplifies the task by allowing the user to specify desired properties directly, with the generation of constraint-satiisfingg expressions entirely managed by the Prolog prover. The implementation of this simplified knowledge base and abduction process detailed the following Prolog code.
X
X Background Knowledge WXXXΕллККииКллМЕЕX
% Declare the dependent variables dep_var ([u, vl).
X Rule for fraction expressions: (X/Y)
X Cur_num Counter for the number of operations used far X 0p_limit: Maximum operations allowed
fraction_ expr (frac Y), Cur_num, 0p_limit) -->
['(],
mitiplication expr (x, Cur _umi, 0p_limit 1), % Parse numerator (X)
[11],
addition_expr(1" Cur Cur_num2, 0p_limit Cur numi 1), % Parse denominator (Y)
(Cur_nue is Cur_numt Cur_numt + 1), % Calculate total operations used (adding for fraction)
())].
X Rule for multiplication expressions: (I+Y) x First clause handles actual multiplication auttiplication.epp (mul (x, Y), Cur_num, 0p-limit) ['],,
unit_expr (x, Cur_numt 0p-limit 1),, % Parse first operand (X)
[+]
multiplication eepeettioneeep(((( (Y, Cur_numt 0p fimii Cur_numi 1), % Parse second operand (Y)
{Cur_num is Cur_numi + Cur_ummm Cur_num2 + 1}, {>+constraint (X, doesn't violate integrity constraints [))].
Y))}, (mul
Check that this multiplication
% Second clause for multiplication.ens
X This is hog the recursion terminates multiplication.eppe
Cur_num, 0p_limit) --> (0 =< 0p_limith,
unit_expr(X, Cur_num, 0p-limit).
X Rule addition expressions: (X+Y) addition_expe (X, Y), Cur_num 0p-limit) -->
[(]]
multiplication (X, Cur_mumt 0p_limit 1), % Parse firat operand (X)
[+],
addition_expr((. Cur_num2, 0p_limit Cur numi 1), %
Parse second operand (Y)
()
{Cur_num is Cur_numi Cur_num2 1).
% Second clause for addition_expi
% This is how the recursion terminates Aldition_expr(e (X, Cur_num, 0p_limit)
20

---

(Stream,, Y), (Stream))), mitippication_exptr (x, Cur_num, 0p_limit).
X Rule for basic units
X This is a terminal rule that matches a single variable from the dependent variable list unit_expr(X, Cur_num, 0p_limit) -->
Cur_num is 0, Cur_num -< 0p_limit,
var (DepVarList), member (x, DepVarList) X x must be one of the dependent variables
}.
[x]
X Abducible Predicates /
X Decide which terms can utilized to explain the observation, i.e., which values of X can satisfy expression_term (X, 0p_limit)
expreasion_teras (X, 0p_limit) expression_texs wermmonttteteeee (x, 0p_limit) --> multiplication.expt (x, 0p.limit).
craction_expr(s 0p_limit).
Integrity Constraints 4XXXIHX4XXXXXXXXIHHG
X Constraint to enforce a canonical ordering of multiplication terms
X This prevents generating both (a+b) and (b+a) which mathematically equivalent
constraint (X, f)) :- extract_firs.mml.cea
FiratTerm),
X 0> FirstTerm
X Helper predicate extract first term in multiplication chain
% For (ul mul(b, c)), it returns 'a'
% For a simple term like 'u', it returns ' (Term, FirstTerm) ;- multiplication
fxtrat_ffrst_mu_teet
functor (Term, mul, 2) % Check if Term is
>> arg Term, FirstTerm) If yes, extract first argument
FirstTerm Term % If not, the term itself is the first
% Utility predicate to write all generated expressions to file write_to_fill open opont'uut.txt: write, Stream),
findall ([x, Y], pbrase(oxpressionn.co termii (X, 2), Y), L), % Generate all expressions with max 2 operators forall (member Y], L), (srite(Stream, I), ni(trream),
close (Stream).
C A Brief Introduction to POSS
POSS leverages Pareto Optimization to solve the following subset selection problem:
21

---

Algorithm 2 POSS
Input: all variables = = (, Xn..
given criterion and an integer parameter [1,
Parameter: the number of iterations 7" and an isolation function / {0, 1)" + IR
Output: subset of V with at most & variables
1: Let {0}m and P {s}.
2: Let = 0.
3: while t < T do
4: Select from uniformly at random.
5: Generate s from by flipping each bit of s with probability 1/n.
6: if 72 E such that T(2) T(s') and
(z.01 &0.01 A 2.02 S 8 8'002) or (z.01 S 58.01 A.02 < <.
then
7:
Q =
E 1()) A2 2.02).
8:
P= P (s')
9:
f=t+1.
10:
Return arg c s f(s)
Definition 6 (Sparse Regression), Given all observation variables V {Xi..... Xn), a predictor variable Z and a positive integer k, define the mean squared error of subset F C V as
(
( (Z-D
27
MSF2,F = min E
x;
aERISI
iEF
Sparse regression is to find a set of at most & variables minimizing the mean squared error, i.c.,
arg min MSEZ,F PCV
s.t.
[F]
We use it to roughly select the subset that contains possible relevant terms. The following is the pseudo-code of POSS:
In Algorithm 2 01 represents residual and 0g represents the subset size. : {0, 1)"h + determines if two solutions are allowed to be compared: they are comparable only if they have the same isolation function value.
D More Details on FEMTrain
D.1 A Brief Introduction to FEM
Finite Element Method (FEM) is a well-established numerical method that provides a systematic framework for approximating solutions to PDEs.
We use scalar field u with one spatial dimension, [0, LJ. as a showcase to explain the process of using FEM to solve the following type PDE:
0uu f (u, x. O 0 ........
(8)
The first step to apply FEM is to discretize the spatial domain into smaller, finite subdomains called elements. Given ordered set of points X = {21 ,,, .n with T₁ 0 and 2n, L, they split the domain into n - 1 elements: [x1, .2l.--. ., [In-1, The second step is to construct set of
basis functions {dh,, .2,,.., Φn) approximate the solution u. di is defined follows:
0,
if < 8i-1 or x > Ti+1.
(x(x))
x-
if fi < x < Ti
if & S Tit1.
0141-27
22

---

The solution u is then approximated by ū defined as:
④ (t,
Σ 0. (1()) (x),
where (2(44 = v(t, w4). The third step to transform Eq. (8) into weak formalation by multiplying both sides with each basis function and integrating over the domain:
Je d.u oi dx = L' (u,,I.... dæ.
(10)
Substitute Eq. ( into Eq. 10 get linear system:
M (1)) F()), dt
(11)
where My = Jo Φ: o, dr is the so-called mass matrix in finite element analysis, at()) is the vector of coefficients in Eq. ( and F.(t) = for (ū, x....) o, dz. Eq. (IT) an Ordinary Differential Equation (ODE) system of a. Therefore, given the initial state, can calculate F and forward a iteratively, thus forecasting ü at any time.
D.2 FEM-based Residual
Firstly, we discretize the spatial domain ! into a set of cubic elements T = {2, C R Let X = (x, E )- denote their vertices, and = {Φ. Ω + R for the corresponding basis functions. Secondly, we use & to approximate 8,u(t, x) & f( (u(t, x), x) and expression terms {f (u(t, x), x )=) identified by MSS as follows:
N
F(0) (u(t, x) = Σ 10 (1) (x), 0, 0,,,, NE,
(12)
where ( (t) f x;), x;)). which can be evaluated through the pre- pre-trained neural network.
Thirdly, we transform (0) ENer & p(0) into the weak formulation:
(f())
Σ Ea
1,2.
(13)
where & is the coefficient as introduced in Eq. ①.. and (u, 0) Sg dæ. Substituting 12 into 13) at a given time point (omitted for brevity) yields the linear system:
Ma (0) MAC.
(14)
where Mij (0h,, Φ3); (1)
(al (9)
2).
(1 (2) a (Nr)) and &
(5o, .1....,En. )T. Since (1) residual can be utilized optimize neural network parameters and PDE coefficients E:
No are outputs of the neural networks, the following HEM-based
(M (9) AE)
(15)
Replacing Res(NNg, F) in Eq. 2 with this residual gives the loss function utilized in FEMTrain.
E Proofs
E.1 Theorem
Lemma 1. Let é1 and ea be two expressions generated by <sub>, Then, e₁ e2 can be rewrillen as an expression generated by <sub>.
23

---

Proof. According to Definition 2 can be expressed as b, e2 can be expressed as d, where a, b, c, are expressions generated by <add>, <mul>, <par>, <dep>, or <indep?. Without loss of generality, we only need to consider the <add> case, as all other cases are simply degenerate forms of this one. Since (a b) x (c d) = (ac+ + bd) (be + ad), we only need to prove ac + bd and be + ad can be rewritten as expressions generated by <add>. Since and are summations of multiplication or unit expressions, it's obvious that we can rewrite ac as an expression generated by <add> according to basic algebraic laws, as bdl bc, and ad. Therefore, (ac + bd) (bc + ad) is subtraction between two expressions generated by <add>, which can be rewritten as an expression generated by <sub>.
Lemma 2. Let and be two expressions generated by <div>. Then, 62 can be rewritten as an expression generated by <div>.
Proof. Without loss of generality, we assume ε = -i e- where and are expressions generated by <mul>; b, c, e, and f are expressions generated by <add>. According to Lemma [ the denominator of 61 e2 be written expression generated by <sub>. The numerator of e1 x 62 can be rewritten as an expression generated by <mul>, thus 61 x e2 can be rewritten as an expression generated by <div>.
Lemma 3. Let 61 and e2 be two expressions generated by <div>. Then, the numerator and denominator of e1 op ez, where op E {+,-,/} can all be rewritten as expressions generated by <sub>.
Proof. Without loss of generality, assume é 62 c-f- where and d are expressions generated by <mul>; b, c, e, and are expressions generated by <add>, We prove the "+" case and others are similar. Since e1 + e2 = =(e-()+d(--)) according to Lemma ale - f) + d(b c) and (b-c)(e-f)
(b (b-c)(ee - f) are all expressions generated by <sub>.
Lemma 4. Let e be an expression generated by <div>. Then, de can be subsumed by a set H' C H. where x an arbitrary independent variable.
Proof. Without loss of generality, we assume bー where is an expression generated by <mul>; and are expressions generated by <add>, According differentiation rule, 0g g (ba (b-) (b-c), According Lemma (b c)² an expression generated by <sub>, According (6 (b-c)²
to the differentiation rule and basic laws of multiplication, the numerator is an expression generated by <sub>. Thus, can be subsumed by a set H' H.
Proof of Theorem 1
We use structural induction prove Theorem |
1. Base Case: For each expression term s with zero operator, there exists a term h E H such that h.
This case is obviously true.
2. Inductive Hypothesis: Assume that for each expression term with no greater than operators, there exists subset of H' C H such that H' subsumes s.
3. Inductive Step: Show that for each expression term with + operators, there exists subset of H' H such that H' subsumes s.
If 8 is an addition or subtraction, according to the inductive hypothesis, the conclusion is obviously true.
If s multiplication, assume s a x b, according to the inductive hypothesis, a and can be subsumed by two sets H,, C H and H, C H respectively. Without loss of generality, we assume each element of H, and H, are all expressions generated by <div>. According to the distributive law of multiplication and Lemma the conclusion holds.
24

---

If s is division, let s = a/b, according to the inductive hypothesis, and b can be subsumed by two sets Ha H and H, H respectively. Without loss of generality, we assume each element of Ha and H, are all expressions generated by <div>. Then, we only need s+d to prove c+] can be subsumed by a set H' C H, where c, d, e, and f are all expressions generated by <div>. Other cases are degenerated eases. According to Lemma e+f can
rewrite as g/h where g. h, j. and h are all expressions generated by <sub>. According 37R
to Lemmal it's obvious that there exists a set H' H subsumes gk
If s is a partial derivative, let s de. According to the inductive hypothesis, e can be subsumed by set H. H. According to the linearity of differentiation, we can assume is single element of Without loss generality, we assume an expression generated by <div>. According to Lemma + the conclusion holds.
4. Conclusion: By the principle of mathematical induction, Theorem holls.
E.2 Proof of Theorem 2
Proof. Let Fi be subset of F with cardinality 2 that realizes the -leell margin, such that LevelMargin (k1) = (FF).
Since b l2 S [F). we construct new set F2g by adding l2 k arbitrary terms from F \ Ff to FY. By construction, F S F32 and FF = k2.
The minimal residual is found by solving a linear least squares problem. The solution space for the coefficients & corresponding to the terms Ff a subspace that for the terms in Pj. Therefore, the minimal residual over the larger set F22 cannot be greater than the minimal residual over its subset Ff:
|Res" (NN,, FT ll2 2
(Res"(NNo. F22 112
From the definition of Margin (Definition 3p. this directly implies that:
M FF) 2 Margin,
Furthermore, by the definition of k-Level Margin (Definition 4), the he-level margin the minimum margin over all subsets of size k2. The margin of our constructed subset F2J must therefore be greater than or equal this minimum:
Margin, (F)) LevelMargin,(xs).
Combining these inequalities, we arrive at the desired result:
LevelMargin = , 2 Margin, 2 LevelMargin (k2))
This completes the proof.
More Experimental Details
F1
Equation Form
Table 3: Summary of Governing Equations
Equation
Form
Burgers' Equation
Ozu = uZ,u + (0.01) :zz
(0,u -0.50 v(u² + v2)
Schrödinger Equation
0₂v = 0.50 u + u(u² + 2²)
(U,u -ud,r u v0, u 0,v' 0.01 (8 + uy
Navicr-Stokes Equation
0,v = -udzv vZ,y 0,p + 0.01 (0zz v + 0yy
25

---

F.2 L2 Relative Error
Definition 7 (L2 Relative Error). Let f [f fe..... ful be vector of true values and vector of approximations. The L2 Relative Error (L2RE) is defined as:
[f1.f.... fe...., fn]
|l flle
VE ff fP?
L2RE
lifll
VE
" i=1 [A
Please note that for coefficient estimation, the 1.2RE simplified as follows:
The benefit of using this metric is that errors in large coefficients do not dominate and obscure errors in small coefficients.
F.3 Network Architecture and Training Details
Our model's architecture is a multilayer perceptron (MLP) composed of hidden layers. For the Burgers' Equation, each layer has 20 neurons. This architecture is adapted for the Schrödinger Equation by modifying the output layer predict variables (the real and imaginary parts of the solution). For the more complex Navier-Stokes Equations, we increase the network's 's width to 40 neurons per hidden layer.
The training data consists of x 10 (40%), x 10a (40%), and 2 x 10 (20%) measurement points for the Burgers', Schrödinger, and Navier-Stokes experiments, respectively. We hold out 20% of the training data as a validation set, which is used to select the derivative computation weights and the hyperparameter A. The networks are pre-trained for 5,000 epochs for the Burgers' and Schrödinger equations, and 100,000 epochs for the Navier-Stokes equations. Across all tasks, we use a consistent initial sparsity of 10 and margin threshold of 0.05. During consistency optimization, MSS is applied every 1,000 epochs to refine the discovered equation. The values for À are selected based on validation error and are summarized in Table All experiments were conducted on a server equipped with four NVIDIA A6000 GPUs, each with 48 GB of memory.
Table 4: Hyperparameter À settings.
Method
Burgers' Equation
Schrödinger Equation
Navier-Stokes Equations
ABL-PDE (w/o FEM) ABL-PDE
0.05
0.5
200
G
Detailed Experimental Results
G.1
Redundant Terms
Table detaails the specific redundant for the margedd (f) results presented in Table
Table 5: Redundant terms of the dagger dagger-marked marked results in Table 2
Method
Equation
Noise Level
Redundant Terms
Ue Maarr U. WU, u.
ABL-PDE (w/o CO)
Schrödinger Eq.
All Levels
0y:
0%
ug:
Ve: Un
DL-PDE++
Navier-Stokes Eq.
ue: p²u, Var
5%
ve: Py puty
ue: p² u, vu p² pppp vc U, puay Ux
10%
26

---

G.2 Test Error
Table presents the detailed results test errors on three benchmarks. The area of the shaded circles in Figure is proportional to the square root of the standard deviation. We use the square root to improve the visibility of the standard deviation of ABL-PDE.
Table 6: Test error comparison results three benchmarks. The evaluation metric is the sum of 1.2RE (%) on each dependent variable. Results of ABL-PDEE (w/o FEM) and ABL-PDE are the mean and standard deviation of experiments with different hyperparameters.
Burgers' Equation
Schröddinger Equation
Naiver Stokes Equation
Method
ABL-PDE (wlo CO)
1,79
1.83
1.92
2.25
2.37
3,93
2.85
3,38
4.70
(wlo
4.02
ABL.PDE
0.68 0.01 0.72 + 0.02 0.85 0.03
1.22 0.01 1.36 + 0.01 1.75 0.01 2.64 0.05 3,17 + 0.04 3.8
G.3
Redundancy Reduction
Table presents the detailed number of candidate terms.
Table 7: Number of candidate terms w.r.t. operator number limit.
Number 2
10
Full
105
44205
7.8E+9
2.4E+20
2.4E+41
2.3E+83
2.1E+167
1.7E+335
1.2E+671
Canonical
39
143
510
1784
6141
20849
69948
234789
764761
ABL-PDE
11
29
69
165
387
879
1955
4254
9053
27


Reviews:
[{"review_id": "hCj0ZXVG37", "reviewer": "Reviewer_DouY", "summary": "The authors present a new approach for automatic discovery of partial differential equations from data using neural network learning and first order logic based reasoning. Current approaches for partial equation discovery suffers from several issues including - large hypothesis space which makes it challenging to search for the correct equation, inaccurate derivative estimation during expression evaluation as only discrete measurements are available. To address these issues, the authors propose to use first-order-logic-based reasoning to constrain the hypothesis space of possible expressions/equations, and neural network based derivative estimation for expression evaluation. The authors’ framework also contains an optimization module to ensure consistency between these two modules. Experimental results are provided to demonstrate the efficacy of the proposed framework.", "questions": "Mentioned in the strength and weakness section", "limitations": "Please see the strength and weakness section", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strength:\nThe authors presented a good introduction and motivation for their work. Related works are discussed along with their limitations. Algorithm for Monotone subset selection is presented in detail. Theoretical analysis is presented. Experimental results are organized with key questions mentioned at the beginning. Ablation study is provided to show the importance of different modules of the framework.\n\nWeakness:\n1. The authors may consider adding the following details to their manuscript:\n2. Why is it necessary to have a consistency maximization process? There can be one/two sentences between line 54-62 describing the intuition behind it.\n3. 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?\n4. 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?\n5. 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.\n6. 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.", "was_revised": true, "final_justification": "I had questions about the manuscript, the authors clarified and provided more details in their response. I have also read other reviewer comments, and I think most of the major concerns have been addressed by the authors. Therefore, I will suggest accept for this paper.", "rebuttal": "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.\n\n---\n\n**Q1.** \"Why is it necessary to have a consistency maximization process?\"\n\n**A1.**\nThanks 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.\n\nWe 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.\nThus, 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.\nIn 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.\n\n---\n\n**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?\"\n\n**A2.** \nThanks 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**.\n\nSpecifically, 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.\n\nThis 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.\n\nReferences:\n\n[1] DLGA-PDE: Discovery of PDEs with incomplete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 2020.\n\n[2] Deep-learning based discovery of partial differential equations in integral form from sparse and noisy data. Journal of Computational Physics, 2021.\n\n[3] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). Physical Review Research, 2022.\n\n[4] Weak SINDy for partial differential equations. Journal of Computational Physics, 2021.\n\n[5] D-CIPHER: discovery of closed-form partial differential equations. NeurIPS, 2023.\n\n---\n\n**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?\"\n\n**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.\n\nWe 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.\n\nExperimental results demonstrate the effectiveness of consistency optimization in improving the accuracy of derivative estimation.\nAs 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**.\n\n---\n\n**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.\"\n\n**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.\n\nFirst, we present a vanilla Prolog program without operator number limit and integrity constraints.\n```prolog\n% Logical facts\nunit_expr(u).\nunit_expr(v).\n\n% Logical rules\nmultiplication_expr(X) :- unit_expr(X).\nmultiplication_expr(mul(X, Y)) :- unit_expr(X), multiplication_expr(Y).\n```\nThe 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`.\nWhen 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.\n1.  **First Rule: Base Cases**\n    * Prolog first tries this rule. It looks for simple \"unit\" expressions.\n    * It finds `unit_expr(u)` and reports `X = u`.\n    * Then it finds `unit_expr(v)` and reports `X = v`. These are the simplest, non-nested expressions.\n\n2.  **Second Rule: Recursive Cases**\n    * After the simple cases, Prolog moves to the second rule to build more complex expressions.\n    * 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`.\n    * **How it works:**\n        * It picks a `unit_expr` for the first part.\n        * 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.\n        * This recursion allows it to build nested structures like `mul(u, u)`, then `mul(u, v)`, `mul(u, mul(u, u))`, and so on.\n        * **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.\n\n**Operator Number Limit**\n\nOne 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.\n\n```prolog\nunit_expr(u).\nunit_expr(v).\n\nmultiplication_expr(X, Op_limit) :- 0 =< Op_limit, unit_expr(X).\nmultiplication_expr(mul(X, Y), Op_limit) :- 1 =< Op_limit, unit_expr(X), \n                                            multiplication_expr(Y, Op_limit - 1).\n```\n\nThe 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.\nWhen querying `?- multiplication_expr(X, 1).`, only `u`, `v`, `mul(u, u)`, and `mul(u, v)`, `mul(v, u)`, `mul(v, v)` will be generated.\n\n**Integrity Constraints**\n\nSince we do not want to both generate `mul(u, v)` and `mul(v, u)`, we introduce the following integrity constraints to the second rule.\n\n```prolog\nmultiplication_expr(mul(u, v), Op_limit) :- 1 =< Op_limit, unit_expr(u), \n                                            multiplication_expr(v, Op_limit - 1),\n                                            \\+ integrity_constraint(mul(X, Y)).\nintegrity_constraint(mul(X, Y)) :- Y = mul(Z, _) -> X @> Z ; X @> Y.\n```\nThe semantics of `\\+` is \"not,\" meaning that for a valid `mul(X, Y)`, the boolean value of `integrity_constraint(mul(X, Y))` should be `false`.  \nThe 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.  \nBy introducing this integrity constraint, we can prevent the generation of terms like `mul(v, u)`.\n\nFor a more comprehensive example of our logic-based knowledge base implementation, please refer to Appendix B."}, {"review_id": "rjRfdPfG7z", "reviewer": "Reviewer_htDg", "summary": "The idea of this paper is to integrate a neural network for fitting the solution field (and estimating derivatives) with a first-order logic (FOL) knowledge base that generates canonical expression terms via abductive logic programming.", "questions": "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.", "limitations": "They test only on well-known equations (Burgers, Schrödinger, Navier-Stokes in 2D).\nIf the initial pre-training fails (very sparse or chaotic data), the abductive loop may not recover.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:  Combining FOL abductive learning to systematically constrain symbolic PDE search is quite original.\nWeaknesses\nWhile the canonical forms reduce combinatorial explosion, the abductive reasoning and MSS steps could still become heavy for large systems. It would help to see timings or scaling on more complex multi-physics PDEs or with >2 dependent variables.\nThough 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.", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer htDg for insightful review and appreciation on the originality of our work. Below, we answer your questions in detail.\n\n---\n\n**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.\"\n\n**A1.**\nThanks 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.\nIt 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.\n\n---\n\n**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.\"\n\n**A2.** \nThanks for your suggestion. We highly agree that testing our method in more complex scenarios would better demonstrate its practical applicability.\nHowever, fitting unknown and complex kinetics often involves learning parameters within nonlinear terms.\nAs we mentioned in A1, learning these parameters significantly increases computational cost and requires more advanced optimization techniques, which we consider as future work. \n\n---\n\n**Q3.** \"If the initial pre-training fails (very sparse or chaotic data), the abductive loop may not recover.\"\n\n**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.\n\nIn 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.\n\nFor 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.\n\nReferences:\n\n[1] DLGA-PDE: Discovery of PDEs with incomplete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 2020.\n\n[2] Deep-learning based discovery of partial differential equations in integral form from sparse and noisy data. Journal of Computational Physics, 2021."}, {"review_id": "LKQdmoOa8T", "reviewer": "Reviewer_dxqN", "summary": "This paper introduces ABL-PDE, a novel abductive learning-based approach for PDE discovery. The key contribution is the transformation of the traditional unidirectional discovery pipeline into a bidirectional enhancement loop, improving both structure and parameter identification. Specifically, the method incorporates FOL KB to represent a broad class of PDEs with reduced redundancy and enhanced interpretability. Extensive experiments across multiple PDE tasks and noise levels demonstrate that ABL-PDE achieves state-of-the-art performance in both equation structure recovery and coefficient estimation.", "questions": "1. 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:\n(1) I noticed that some L2RE values are exactly the same across different methods for the same PDE and noise level. Is this because the same Ridge regression is used in the final step of those methods?\n(2) Why are standard deviations (± values) reported for some results but not others? Is this because those experiments are deterministic, with no variation across runs?\n(3) I understand that the dagger symbol indicates the discovered PDE includes redundant terms. This is common in symbolic regression. However, how severe are these cases? For example, one extra term with a tiny coefficient might be fine, but having five extra terms with meaningful coefficients in a three-term PDE could be considered a failure. Also, does the dagger indicate that all correct terms are present and only redundant ones are added, with nothing missing?\n(4) If an extra term is included in the discovered PDE, how does the L2RE be calculated? Do you ignore the coefficients of redundant terms and only use the correct terms when computing the error?\n\n2. Could the authors compare their method and its variants with more recent and representative baselines? PDERidge only applies Ridge regression on a known PDE structure, which is not a complete PDE discovery method. DL-PDE++ (2021) is somewhat outdated, though still acceptable. 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).\n\n3. Since the overall optimization process alternates between MSS and FEMTrain, can the authors provide a proof or argument that this process always converges?\n\n4. The caption of Figure 3 mentions that 5 different values of lambda (used in Equation 2) are considered. However, neither Figure 3 nor Appendix G clearly shows or discusses these different lambda values. Could the authors clarify this?\n\n5. In section 6.4, it is stated that the lambda values used are [0.5λ, 0.75λ, λ, 1.25λ, 1.5λ]. Does this mean five separate configurations with λ set to 0.5, 0.75, 1.0, 1.25, and 1.5 respectively? The wording is a bit unclear.", "limitations": "Yes.", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strengths\n1. The authors replaces the traditional unidirectional pipeline with a feedback loop between symbolic and numerical components, enhancing robustness and accuracy.\n2. The authors introduce a knowledge base to constrain the candidate search space, effectively reducing the problem size to a computationally manageable scale. This approach is both practical and innovative - I like this idea.\n\n### Weaknesses\n1. The Experiments section includes a limited number of baseline methods for comparison. The method could benefit from comparisons with more recent or diverse PDE discovery methods.\n2. The alternating optimization between MSS and FEMTrain lacks a theoretical guarantee of convergence, which may affect reproducibility or stability in some scenarios.\n3. The choice of using L2RE of term coefficients can not clearly reflect the performance of equation discovery.", "was_revised": true, "final_justification": "I have read the author rebuttal and the discussions among all reviewers. The authors have addressed my concern regarding hyper-parameter settings. However, my core concern remains - the choice of reporting only the L2RE of term coefficients as the sole experimental metric is insufficient. I recommend also reporting the recovery rate (success rate) together with the equation value MSE (or RMSE, or any MSE-based metric) calculated on the entire equation, not just the coefficients.\n\nThe quality of equation discovery depends both on the proportion of correct terms recovered and on the accuracy of the predicted values. This is also why the reported performance improvement of ABL-PDE over the baselines remains unclear to me (e.g., baseline performance comparison in A2), as the chosen metric is not well-suited for a comprehensive evaluation. For the baseline experiments mentioned in all reviewers' comments, if available, please include results using the updated metrics in the revised version.\n\nFurthermore, as noted in Q3 and A3, a convergence proof between MSS and FEMTrain is necessary to establish a solid foundation for ABL-PDE. Taking these points into account, I will maintain my rating for this paper.", "rebuttal": "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.\n\n---\n\n**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?\"\n\n**A1.** \nThanks for your helpful question. We will revise our paper to clarify these points.\n1. Exactly. We used the same pre-trained neural network and the same regression method to ensure a fair comparison.\n2. 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.\n3. 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):\n\n|Method|Burgers' Equation|Schrödinger Equation|||Navier-Stokes Equations|||\n|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|\n||10%|0|5%|10%|0|5%|10%|\n|DL-PDE++|||||1|4|7|\n|ABL-PDE (w/o CO)|2|3|3|3||||\n\n4. Exactly. When computing the error, we ignore the coefficients of redundant terms and only use the correct terms.\n\n---\n\n**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).\"\n\n**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.\nSince 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.\nSimilarly, 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.\nWe 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.\nDue 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}$.\n\n|  | Burgers' Equation | KdV Equation | Chafee-Infante Equation |\n|---|---|---|---|\n| Ground Truth | $u_{t} = −uu_{x} + 0.1 u_{xx}$ | $u_{t}=-0.0025u_{xxx}-uu_{x}$ | $u_{t}=u_{xx}-u+u^3$ |\n| SGA-PDE[1] | $u_{t} = −1.0011 uu_{x} + 0.1024 u_{xx}$ | $u_{t}=-0.0025u_{xxx}-1.0004uu_{x}$ | $u_{t}=1.0002u_{xx}-1.0008u+1.0004u^3$ |\n| ABL-PDE (w/o CO) | $u_{t} = -0.83uu_{x}+0.096u_{xx}-0.026u+0.038 u^3$ | $u_{t}=-0.00237u_{xxx}-0.988uu_{x}$ | $u_{t}=0.843u_{xx}-0.852u+0.916u^3-0.0619uu_{x}+0.182u_{x}$ |\n| ABL-PDE | $u_{t} = −0.999uu_{x} + 0.10002 u_{xx}$ | $u_{t}=-0.00249u_{xxx}-0.995uu_{x}$ | $u_{t}=0.994u_{xx}-0.987u+0.996u^3$ |\n\n|  | Lorenz System |\n|---|---|\n| Ground Truth | $\\dot{x} = −10x+10y$ &nbsp; $\\dot{y} = 28x-y-xz$ &nbsp; $\\dot{z} = -\\frac{8}{3}z+xy$ |\n| SPL[2] | $\\dot{x} = -9.966x + 9.964y$ &nbsp; $\\dot{y} = 27.764x - 0.942y - 0.994xz$ &nbsp; $\\dot{z} = -2.655z + 0.996xy$ |\n| ABL-PDE (w/o CO) | $\\dot{x} = -9.807x+9.956y+0.0323xy-0.0185yz$ &nbsp; $\\dot{y} = 27.72x - 0.976y - 0.989xz$ &nbsp; $\\dot{z} = -2.698z + 0.950xy+0.0238yz$ |\n| ABL-PDE | $\\dot{x} = -10.003x+10.002y$ &nbsp; $\\dot{y} = 27.869x - 0.976y - 0.995xz$ &nbsp; $\\dot{z} = -2.667z + 0.999xy$ |\n\nReferences:\n\n[1] Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). Physical Review Research, 2022.\n\n[2] Symbolic Physics Learner: Discovering governing equations via Monte Carlo tree search. ICLR, 2023.\n\n---\n\n**Q3.** \"Since the overall optimization process alternates between MSS and FEMTrain, can the authors provide a proof or argument that this process always converges?\"\n\n**A3.**\nThanks 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.\n\n---\n\n**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.\"\n\n**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:\n\n|  | Burgers' Equation | Schrödinger Equation | Navier-Stokes Equations |\n|---|---|---|---|\n| ABL-PDE (w/o FEM) | 0.05 | 0.05 | 1 |\n| ABL-PDE | 5 | 0.05 | 200 |"}, {"review_id": "XlRjd1KgfE", "reviewer": "Reviewer_B51F", "summary": "This paper presents a methodology for discovering the differential operator of Partial Differential Equations using Abductive Learning. The method is structured in two steps, first a neural network is used to approximate the given data such that there is no need for an integrator and second a process that contains an outer and an inner loops is considered. The outer loop increases the length of the differential operator, the inner loop picks the candidate operators, use Mean Subset Selection to find the shorted expression and fine-tune the neural network. The process is performed until convergence on a level margin criterion. The method is tested on three PDEs for different levels of noise against different baselines, and on two PDEs specifically against D-CIPHER which cannot handle systems.", "questions": "- Could you provide more information about how the method handles additive coefficients?\n\n- When performing re-writes of complex terms in the grammar, is there a way to derive a bound to this growth such the method remains efficient when considering more complex expressions (ones with forcing terms for example).\n\n-  Maybe I missed it but could you clarify if the candidate terms are the same for all PDEs or not? And if not how they are chosen?\n\n- How many times did you run D-CIPHER and what was the best result?\n\n- Why is the discovery more accurate for 5% error?", "limitations": "There is only one limitation stated in the conclusions. I believe that the authors should include bullets from the weaknesses section to the limitations.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\n- This is an interesting method that combines canonical grammars that are expressive enough, see Theorem 1, together with a consistency optimization where is showed that the k-level margin decreases monotonically. The method essentially proposes the idea of reducing the hypothesis using a first-order logic knowledge base and integrity constraints, and delivers state-of-the-art accuracy for PDEs with different levels of noise. \n\nOverall, it provides a novel method to cut down the hypothesis space of combinatorial optimization methods, which is very valuable. Also the paper is mostly well structured, well written, and the authors provide helpful information in the appendix. \n\nWeaknesses:\n\nClarity: I believe that there are some parts that are not completely clear to me. For example:\n\n- How are constants treated in the Grammar. The non-terminal <unit> allows for (in-)dependent variables, and partials only, meaning that constants, for example a source term \"+2\", cannot be generated. However, constants appear in the target equations, for example 0.01/pi in the Burger's. In my understanding, the authors rely on the \\xi vector in the FEM training to account for that but this doesn't work for additive constants, such as u_t - D u_xx - S where S a constant source. \n\n- There is the following subtle issue in the method (at least to my understanding): The strength of using a grammar is to keep the library of candidate term small such that the MSS choice remains tractable. Theorem 1 states that any expression, even if it is nested, can be re-written as linear combination of terms from the canonical grammar. However, when this happens a single complex node is substituted by multiple simpler nodes, which means that you might lose the benefit of restricting the space of candidate expressions and having a large library again. \n\n- I believe that there is one more subtle issue with Theorem 1. In my understanding, Theorem 1 guarantees equivalence of symbolic expressions, but does not guarantee that those equivalences hold as PDE operators acting on a solution space with the correct boundary conditions. The difference is that without discussing boundary condition assumptions, probably smoothness also because maybe the solution has weak derivatives, you can be sure that the expansion holds in the weak form where you consider in the FEM residual. \n\n- In my understanding, Theorem 1 assumes that derivatives and variables commute and the grammar contains the basic terms, but it seems like constants are missing from the grammar. \n\n- It is not clear how the FEM based training handles boundary conditions in the weak formulation. It seems like the training ignores what happens on the boundary. Moreover, I believe that there should be some discussion on how the conditioning of the mass matrix affects the overall  training, as M might be ill-conditioned for coarse meshes. \n\n- It is not clear to me how expensive and numerically stable is the LevelMargin step as it needs to compute orthogonal projections which is expensive and might be numerically unstable for colinear terms.\n\nLemma 5:\n\n- In the proof it is stated that \"Without a loss of generality we assume $v \\perp Y$.\". This statement is not justified and the derivations that follow depend critically on that orthogonality. This statement is not without a loss of generality because if v has a component inside Y, replacing by its orthogonal component v_{ortho} = v - Proj_Y v changes the left-hand side of the inequality because the projection onto span(v) does not coincide with the projection on span(v_{ortho}).  This means that the inequality in Theorem 2 does not necessarily hold\n\n. If the authors impose the condition $v \\perp Y$, the inequality is true but the method does not apply to arbitrary grammars, because you would need to guarantee the orthogonality. \n\n- The results for D-CIPHER are (in my understanding) in contradiction with the original paper, that shows that D-CIPHER is robust to noise.", "was_revised": true, "final_justification": "The reviewers answered my questions and I am satisfied with the responses.", "rebuttal": "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.\n\n---\n\n**Q1.** Clarification on constants in the grammar and additive coefficients.\n\n**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.\n\n---\n\n**Q2.** Concern regarding the re-write of complex terms.\n\n**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.\n\nFirst, 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.\nMoreover, 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.\n\n---\n\n**Q3.** Clarification on the FEM-based residual and operator equivalence in Theorem 1.\n\n**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.\n\nTo 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. \n\nRegarding 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.\n\nThank 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.\n\n---\n\n**Q4.** Clarification on the computational cost and numerical stability of the LevelMargin step.\n\n**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.\nAs 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.\n\n---\n\n**Q5.** Clarification on the proof of Lemma 5.\n\n**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.\n\n**Proof:**\n\nLet $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)$.\n\nSince $\\text{span}(Y \\cup \\{v\\}) = \\text{span}(Y \\cup \\{v\\_{Y^\\perp}\\})$, the RHS remains unchanged:\n$$\n\\\\|\\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.\n$$\n\nFor the LHS, consider the relationship between $v_Y$ and $\\text{span}(X)$:\n\n1. **Case 1:** $v_Y \\in \\text{span}(X)$  \n   If $v_Y$ lies in $\\text{span}(X)$, then $\\text{span}(X \\cup \\{v\\}) = \\text{span}(X \\cup \\{v_{Y^\\perp}\\})$. The inequality becomes:\n   $$\n   \\\\|\\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.\n   $$\n2. **Case 2:** $v_Y \\notin \\text{span}(X)$  \n   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:\n   $$\n   v = v\\_{Y^\\perp} + v\\_{Y,X} + v\\_{Y,X^\\perp}.\n   $$\n   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:\n   $$\n   \\\\|\\text{Proj}\\_{\\text{span}(X \\cup \\{v\\})}(u)\\\\|\\_2 \\ge \\\\|\\text{Proj}\\_{\\text{span}(X \\cup \\{v\\_{Y^\\perp}\\})}(u)\\\\|\\_2.\n   $$\n\nFrom Lemma 5:\n$$\n\\\\|\\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.\n$$\n\nCombining results, we have:\n$$\n\\\\|\\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.\n$$\n\nThis confirms the inequality in all cases.\n\n---\n\n**Q6.** Clarification on the comparison with D-CIPHER.\n\n**A6.** Thanks for your question. D-CIPHER is a SOTA PDE discovery method.\nOne 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.\n\n||Burgers' Equation|||Heat Equation|||\n|---|---|---|---|---|---|---|\n||0|5%|10%|0|5%|10%|\n|D-CIPHER|1.61e-2|8.13e-2|1.49e-1|4.37e-4|1.87e-2|2.62e-2|\n|ABL-PDE|**1.77e-3**|**4.17e-3**|**7.11e-3**|**4.04e-4**|**5.03e-5**|**1.17e-3**|\n\nThe 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.\n\n---\n\n**Q7.** Clarification on the candidate term selection.\n\n**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.\n\n---\n\n**Q8.** Clarification on the discovery result under 5% noise.\n\n**A8.** \nThank you for your detailed review of our paper! This non-monotonic result likely arises from our consistency optimization (CO) stage for two reasons:\n1.  **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.\n2.  **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.\n\nNotably, the physical field fitting error shown in our paper's Figure 3 remains monotonic with noise, which supports the overall validity of our experiments."}]

Author rebuttal:

Meta-review: This submission combines deep learning, a first-order logic knowledge base and consistency optimization to discover symbolic PDEs. The proposed bidirectional loop appears to be novel and the experimental results across benchmarks and noise levels are quite strong. The authors provided an extensive rebuttal and engaged extensively in the discussion throughout which they managed to mostly clarify the theoretical soundness and convinced most reviewers of the originality and applicability of the work. One reviewer remained slightly unconvinced but does not strongly oppose acceptance. Given the convincing clarifications and broad agreement among reviewers, I recommend accept as a poster.
Decision: Accept (poster)