Data-Driven Discovery of PDEs via the Adjoint Method
Authors: Mohsen Sadr
Tony Tohme
KAMAL YOUCEF-TOUMI
Venue: Accepted by TMLR ()

Abstract: In this work, we present an adjoint-based method for discovering the underlying governing partial differential equations (PDEs) given data. The idea is to consider a parameterized PDE in a general form and formulate a PDE-constrained optimization problem aimed at minimizing the error of the PDE solution from data. Using variational calculus, we obtain an evolution equation for the Lagrange multipliers (adjoint equations), allowing us to compute the gradient of the objective function with respect to the parameters of PDEs given data in a straightforward manner. In particular, we consider a family of temporal parameterized PDEs that encompass linear, nonlinear, and spatial derivative candidate terms, and elegantly derive the corresponding adjoint equations. We show the efficacy of the proposed approach in identifying the form of the PDE up to machine accuracy, enabling the accurate discovery of PDEs from data. We also compare its performance with the famous PDE Functional Identification of Nonlinear Dynamics method known as PDE-FIND \cite{rudy2017data} among others, on both smooth and noisy data sets. Even though the proposed adjoint method relies on forward/backward solvers, it outperforms PDE-FIND in the limit of large data sets thanks to the analytic expressions for gradients of the cost function with respect to each PDE parameter.

Published in Transactions on Machine Learning Research (09/2025)
Data-Driven Discovery of PDEs via the Adjoint Method
Mohsen Sadr"
msadr@mil. edu
Massachusetts Institute of Technology, USA
Paul Scherrer Institute, Switzerland
Tony Tohme"
tohme@mit.cm edu
Massachusetts Institute Technology, USA
Kamal Youcef- Toumi
youcef@mit. edu
Massachusetts Institute of Technology, USA
Reviewed OpenReview: https: oppprevve/perrveeeu. net/ forum? id= An3mJ4d1el
Abstract
In this work, we present an adjoint-based method for discovering the underlying governing partial differential equations (PDEs) given data. The idea is to consider a parameterized PDE in a general form and formulate PDE-constrained optimization problem aimed at minimizing the error of the PDE solution from data. Using variational calculus, we obtain an evolution equation for the Lagrange multipliers (adjoint equations), allowing us to compute the gradient of the objective function with respect to the parameters of PDEs given data in straightforward manner. In particular, we consider family of temporal parameterized PDEs that encompass linear, nonlinear, and spatial derivative candidate lerms, and elegantly derive the corresponding adjoint equations. We show the efficacy of the proposed approach identifying the form of the PDE up to machine accuracy, enabling the accurate discovery of PDEs from data. We also compare its performance with the famous PDE Functional Identification of Nonlinear Dynamics method known PDE-FIND Rudy et al. (2017) among others, on both smooth and noisy data sets. Even though the proposed adjoint method relies on forward/backwea /bcckward solvers, it outperforms EDEEFIND in the limit of large data sets thanks to the analytic expressions for gradients of the cost. function with respect to each PDE parameter. Our implementation is publicly available on GitHub
1 Introduction
A significant body of literature on data-driven modeling of physical processes focuses on using neural networks to achieve rapid prodictions from training datasots. The data-driven estimation methods include Physics-Informed Neural Networks Raissi et (2019),
neural networks Eidnes and Lye (2024). , structure preserving Matsubara et al. (2020): Sawant al. (2023), and reduced order modelling Duan and Hesthaven (2024). These methods often provide efficient and somewhat "accurate" predictions when tested as an interpolation method in the space of input or boundary parameters. Such fast estimators are beneficial when many predictions of dynamic system is needed, for example in the shape optimization task in fluid dynamics.
However, the data-driven estimators often fail provide accurate solution the dynamical sys- tem when tested outside the training space, i.c. for extrapolation. Furthermore, given the regression-based nature of these predictors, often they do not offer any error estimator in prediction. Since a wide range of numerical methods already exist for solving classical governing equations, compelling alternative is to
Equal contribution. httpp::/ggtthub.ccmm/m http:://hhhhhccccoocccaaaasaallp

---

Published in Transactions on Machine Learning Research (09/2025)
infer the underlying equations directly from data. Once the governing equation is found, one can either use the standard numerical methods for prediction, or train PINN-like surrogate model for fast evaluation. This way we guarantee the consistency with observed data, estimator for the numerical approximation, and interoperability. Hence, learning the underlying physics given data has motivated a new branch in the scientific machine learning for discovering the mathematical expression as the governing equation given data.
The wide literature of data-driven discovery of dynamical systems includes equation-free modelling Kevrekidis
et al. (2003), artificial neural networks González- Garcia et al. (1998), nonlinear regression Voss et al 1999),
empirical dynamic modeling (Sugihara et al. (2012); :
Ye et al. (2015). modeling emergent behavior Roberts
(2014), automated inference of dynamics Schmidt et al. (2011 Daniels and Nemenman (2015a bb. normal
form identification in climate Majda et al.
(2009) nonlinear Laplacian spectral analysis Giannakis and Majda
(012)). modeling plasma physics Alves and Fiuza (2022). and Koopman analysis Mezić (2013) among others.
There has been significant advancement this field by combining symbolic regression with the evolutionary
algorithms Bongard and Lipson (2007): Schmidt and Lipson (2009); Tohme et al. (2022). which enable the
direct extraction of nonlinear dynamical system information from data. Furthermore, the concept of sparsity
Tibshirani (1996) has recently been employed to efficiently and robustly deduce the underlying principles of
dynamical systems Brunton et al. (2016); Mangan et al. (016).
Related work. Next, we review several relevant works that have shaped the current landscape of discovering PDEs from data:
PDE-FIND Rudy et al. (017)) This method has been developed to discover underlying partial differential equation by minimizing the Lş-norm point-wise error of the parameterized forward model from the data using sparse regression. Estimating all the possible derivatives using Finite Difference, PDE-FIND constructs a dictionary of possible terms and finds the underlying PDE by performing a sparse search using ridge regression problem with hard thresholding, also known as STRidge optimization method. Several further developments in the literature has been carried out based this idea Champion et al. (2019) ; Kaheman et al. (2020) including Weak-SINDy Messenger and Bortz (202T). In these methods, the size (or dimension) of the data set increases, the PDE discovery optimization problem based on point-wise error becomes extremely expensive, forcing the user to arbitrarily reduce the size of data by resampling, or compressing the data using proper orthogonal decomposition. Needles to say, in case of non-linear dynamics, such truncation of data can introduce bias in prediction leading to finding a wrong PDF.
PDE-Net, PINN-SR, and PDE-LEARN. One of the issues with the PDE-FIND is the use of Finite Difference in estimating the derivatives. This has motivated the idea of combining the PDE discovery task with the PDE estimator that avoids the use of Finite Difference. In methods such as PDE-Net [Lng ct al. 20118 2019), PDE-LEARN Stephany and Earls (2024) or PINN-SR Chen et al. 22021) the scarch for weights/ biases of a complicated neural network as a differentiable data-driven PDE est estimator is combined with the sparse search in the space of possible terms to find the coefficients the underlying PDE. Combining
PDE discovery with data-driven estimation of PDEs makes these mothods more expensive than PDE-FIND
in practice.
Hidden Physics Models Raissi and Karniadakis (0118). This method assumes that the relevant terms of the governing PDE are already identified and finds its unknown parameters using Gaussian process regression (GPR). While GPR is an accurate interpolator which offers an estimate for the uncertainty in prediction, its training scales poorly with the size of the training data set as it requires inversion of the covariance matrix.
Contributions. In this paper, we introduce a novel approach for discovering PDEs from data based on the well-known adjoint method, i.c. PDE-constrained optimization method. The idea is to formulate the objective (or cost) functional such that the estimate function f minimizes the Lžznorm error from the data points f* with the constraint that is the solution Here, we consider a parameterized PDE in general form and the task is to find all the parameters including irrelevant ones. By finding the variational extremum of the cost functional with respect to the function f, we obtain a backward-in-limee evolution equation for the Lagrange multipliers (adjoint equations). Next, we solve the forward parameterized PDE as well as the adjoint equations numerically. Having found estimates of the Lagrange multipliers and solution to the forward model f, we can numerically compute the gradient.
parameterized PDE using the method of Lagrange multipliers.
2

---

Published in Transactions on Machine Learning Research (09/2025)
of the objective function with respect to the parameters of PDEs given data in a straightforward manner. In particular, for a lamily of parameterized and nonlinear PDFs, we show how the corresponding adjoint equations can be elegantly derived. We note that the adjoint method has been successfully used before as an efficient method for uncertainty quantification Flath et al. 2011); shape optimization and sensitivity analysis method in fluid mechanics Hughes et al. (1998); Jameson (2003;; Caflisch et al. (2021) and plasma physics Antonsen et al. (2019); Geraldini et al. (2021). Unlike the usual use of PDE-constrained adjoint optimization where the governing equation is known, in this paper we are interested in finding the form along with the coefficients of the PDE given data.
The remainder of the paper is organized as follows. First in Section 2 we introduce and derive the proposed adjoint-based method of finding the underlying system of PDEs given data. Next in Section 3 we present our results on a wide variety of PDEs and compare the solution with the celebrated FDE-FIND in terms of error and computational /training time. In Section we discuss the limitations for the current version of our approach and provide concluding remarks in Section
2 Adjoint method for finding PDEs
In this section, we introduce the problem and derive the proposed adjoint method for finding governing equations given data.
Assume we are given a data set on a spatial/temporal grid g = - go) with g0) ((x(k) z(5)) k == 1,.., N for the vector of functions f* where h is the spatial index and j the time index with ((N.) T being the final time. Here, (kk) E C R" is spatial position inside the solution domain Ω {lt) denotes the j-th time that data is available, and output is discrete map f 9 + R., The goal is to find the governing equations that accurately estimates f at all points on G. In order to achieve this goal, we formulate the problem using the mothod of Lagrange multiplicrs.
Problem setup.
Adjoint method. For simplicity, let us first consider only the time interval t [t0), ,(++1). Consider a general forward model C[-] that evolves N-dimensional vector of sufficiently differentiable functions f(x, t 160) in i E (601)) and æ E Z where the i-th PDE is given by
c,ffl = 0.f + Σ aad V V [P
(1)
d.p
for i = 1,....,... N, resulting in a system of N-PDEs, i.c. the i-th PDE Li predicts fi- Here, x: [mx1; 22, ......... Fnn] is an n-dimensional (spatial) input vector, and f = [fi, .,,.....ff........ SN] is an N-dimensional vector of functions. We use the shorthand fi fi( t) and f = (xx, l). Furthermore, = [PI,...,PN]] and d d= [di, d2...., dn) are non-negative index vectors such that f2 f(" f2 ... fR where p. for i = I...., N denotes the power of fa and
[P]] (d))
(id [SP 53" (R*)
(2)
is an iterated differential operator acting on fo where ((d) for 1, n indicates d,-th derivative in zj dimension, and a.fi denotes the time derivative of the i-th function. We denote the vector of unknown parameters by a = (tidd.p((. where D represents the domain of all valid combinations of i, d, and p.
Having written the forward model [ as general as possible, the goal is to find the parameters a such that f approximates the data points of f at t= +0i+1) +=((++) given the solution f = f at t t= t(o) To this end, we formulate a semi-discrete objective (or cost.) functional that minimizes the n-o-norm error between what the model predicts and the data f on gli+1), with the constraint that f solves the forward model in Eq. Q i.e.
N
cl = Σ Σ (1) (60)) - -1(00,00))) (0+1))² + | (x ())dd)))))))))) )
+ co 00]1
(3)
where |1-112 denotes Lo-norm, and 6 is the regularization factor. We note that PDE discovery task ill-posod since the underlying PDE is not unique and the regularization term helps us find the PDE with the least possible coefficients.
3

---

Published in Transactions on Machine Learning Research (09/2025)
Clearly, given estimates of f and Lagrange multipliers À (A1, A2 ................................ ,N), the gradient of the cost function with respect to model parameters can be simply computed via
ac (-1) /
, ( (iildadl + ddp
(4)
d.s.d.p
where i 1,... N and d = d1 +.. + dn, where |-| denotes Lı-norm. Here, used integration by parts and imposed the condition that À - 0 on the boundaries of Ω al all time L E [0, T.. The compact support of À is motivated by the fact that we consider boundary conditions as known. In case boundary conditions are parameterized, we need to add another constraint to find its paramcters.
The analytical expression can be used for finding the parameters of PDE using in the gradient descent. method with update rule
OC
Oi,d.p Oi,d.p 1I do
(5)
for i = I,..., NN, where y= B3 min(Ar)) dm) is the learning rate which includes free parameter / and scaling coefficient for each term of the PDE, and dmax = max( for all considered d. Let us also define Pmax = max((pl) as the highest order in the forward PDE model. We note that since the terms of the PDFs may have different scaling, the step size for the corresponding coefficient must be adjusted accordingly. Based on our simulation studies, the gradient of the cost function most sensitive to the highest order terms of the PDE. In Appendix D we give a justification for our choice of the learning rale r).
However, before we use Eq. 4 and ()) need to find X, hence the adjoint equation. This can be achieved by finding the functional extremum of the cost functional C with respect to f. First, we note that the
total variation of C. can be derived as
SC = Σ Σ
, (0tD))8(.00) (611) - Σ (z.) - ((()(())))++))))) (6++1)))5 (000)
(-1))dd
di. v d (al)
+
Ot
d.p
where of: denotes variation with respect to f((x,) t) in E (+5), pli+1)), (t)),)(+11)) and V0,,0,, variation with respect to [a = (K) t= ((++1). In this derivation, we discretized the last integral resulting from integration by parts in time using the same mesh as the one of data g(++D, Here again, we used integration by parts and imposed the condition that X + 0 on the boundaries of Ω at all time i E (100) ((++1) for i = N. Note that f( t) is the output of i-th PDE.
Next, we find the optimums of C (and hence the adjoint equations) by taking the variational derivatives with respect to Ss and ( 20+1),
(104),, i.e.
8C = 0
0%:
f
=
Øt
=
ai... a V (
-
(6)
d.p
and
8C
= 0 = (k)) ((01)))))))))))) = 2( (x)(), ((()))1))))) - (k) (,0++)))))))
(7)
)5 o0
for i= I,..., N and j 0, ...,..,, N₁ 1. We note that the adjoint equation 6 for the system of PDEs is backward in time with the final condition at the time t +0++1 given by Eq. 7 In Appendices A and B provide a detailed derivation of adjoint equation and its gradient. In order to make the notation clear, we also present examples for adjoint equations in Appendix G The adjoint equation is in the continuous form, while the final condition is on the discrete points, i.e. on the grid g(++), In order to obtain the Lagrange multipliers in E [((++), {(6), a numerical method appropriate for the forward and adjoint equation should be deployed. Furthermore, the adjoint equation should have the same or coarser spatial discretization as g(++1) to enforee the final condition
4

---

Published in Transactions on Machine Learning Research (09/2025)
Training with smooth data set. The training procedure follows the standard gradient descent method. We start by taking initial guess for parameters e.g- here we take a = 0 initially. For each time interval first we solve the forward model numerically to estimate f(all), (66t1) f(al))(((1)) given the imitial
&E E [t0),
(6+1)), [(0), condition
f(x ((0) = f (
(8)
Then, the adjoint Eq. 6 is solved backwards in time with the final time condition 7 Finally, the estimate for parameters of the model is updated using Eq. 5 We repeat this for all time intervals j = 0,..., N₂ - until
convergence.
In order improve the search for coefficients and enforce the PDE identification, we also deploy thresholding Blumensath and Davics (2009), i.c. set dg,dpp = 0 if b.d threshold, during and at the end of training, respectively. Thresholding improves the regression by reducing PDE space into smaller intermediate PDEs by dropping out irrelevant terms during training. At the end of the process, the terms whose coefficients survive thresholding form the identified PDE.
where a is uscr-defined
In Algorithm 1 we present a pseudocode for finding the parameters of the system of PDEs using the Adjoint method (a flowchart is also shown in Fig. of Appendix E.. For the introduced hyperparameters, note that 3 in the learning rate needs to be small enough to avoid unstable intermediate guessed PDEs, 6o must be large enough to ensure uniqueness in cases where more than one solution may exist, and the thresholding should be applied only when solution to the optimization is not improving anymore up to a user-defined lolerance "thr- For suggested default values, please see the description of the algorithm. For experimental investigation on impact these parameters for a few examples, see Appendix F
We note that the type of guessed PDE may change during the training, which adds numerical complexity to the optimization and motivates the of appropriate solver for each type of guessed PDE, e.g. Finite Volume method for hyperbolic and Finite Element method for Elliptic PDEs. For simplicity, in this work we use the second-order Finite Difference method across the board to estimate the spatial and Euler for the time derivative with small enough time step sizes in solving the forwardbacckward backward equations to avoid blow-ups due to possible instabilities. See appendix C for the analysis on the numerical error for the estimated adjoint gradient. We note that the adjoint method is most effective when there some prior knowledge of the underlying PDE type, and suitable numerical method is deployed.
Algorithm 1 Finding system of PDEs using Adjoint method. Default threshold a 10-3 applied after Naa = 100 iterations, with tolerances 10-9 and Yhhr 10-6, and regularization factor = 10-12. Input: data f*, learning rate n, tolerance 7, threshold applied after Netr, and co-
Initialize the parameters repeat
for j = 0....,
do Estimate f in E (16), (6+1)] by solving forward model I given initial condition
Solve for in / E (16),20+1) given the backward adjoint Eq. (6) and final condition Eq. from data. Compute the gradient using Eq. Update parameters using Eq. end for
if Epochs > Nehr or convergence in a with Thh then Thresholding: set for all such that Jaul < o end until Convergence a with tolerance y Output:
5

---

Published in Transactions on Machine Learning Research (09/2025)
Training with noisy data set. Often the data set comes with some noise. There are several pre- pre-processing steps that can be done to reduce the noise at the expense of introducing bias, for example removing high frequencies using Fast Fourier Transform or removing small singular values from data set using Singular Value Decomposition.
However, we can also reduce the sensitivity of the training algorithm to the noise by averaging the gradients before updating the parameters. Assuming that the noise is martingale, the Monte Carlo averaging gives the unbiased estimator for the expected value of the gradient over all the data set. We adapt the training procedure by averaging gradients over all available data points and then updating the parameters (see Algorithm 2 and the flowchart in Fig. 18 of Appendix D for more details). Clearly, this will make the algorithm more robust al higher cost since the update happens only after seeing all the data.
Algorithm 2 Finding system of PDEs using the Adjoint method with averaging for the computation of gradients over the data set. Default threshold 10- 3 applied after Nihr 100 iterations, with tolerances y 10 and The 10 6, and regularization factor én 10-12
Input: data f", learning rate 7, tolerance 7, threshold or applied after Nehe: and co-
Initialize the parameters repeat
for 0,..... N₁ N2 do
Estimate f in / E (16), (6+11] by solving forward model 1 given initial condition 8
Solve for in & E [16), (6+1)) given the backward adjoint Eq. (6) and final condition Eq. Compute the gradient g0 = =cC /0ax using Eq. 44)
7 from data.
end for
Average the gradient R[aC/Oa] = E, g00) /N2
Update paramcters using Eq. (5) and E(ac/Oa]
if Epochs > Nehr or convergence in with Tehr then Thresholding: set for all such that Jasl < end
until Convergence in a with tolerance 7 Output:
3 Results
We demonstrate the validity of our proposed adjoint- based method in discovering PDEs given measurements on a spatial-temporal grid. We have compared our approach to PDE-FIND in terms of error and time to convergence. All the results are obtained using a single core-thread a 2.3 GHz Quad-Core Intel Core i7 CPU. In this paper, we report the execution time т obtained with averaging over independent runs and use error bars show the standard deviation of the expected time, i.e. derror-bar VB[(T (B((T-B7)) E.)]] Furthermore, also report the true positivity ratio Lagergren al. (0200 defined by
TP
TPR = TP+EN+FP
99)
where TP denotes the number of correctly identified non-zero coefficients, FN the number of coefficients falsely identified as zero, and FP is the number of coefficients falsely identified as non-zero. Identification of the true PDE form results in TPR 1.
3.1 Considered PDEs
In this section, consider the PDE discovery task given data from numerical solutions to variety of problems, including the Heat, Burgers', Kuramoto Sivashinsky, Random Walk, and Reaction Diffusion equations, summarized in Table [
6

---

Published in Transactions on Machine Learning Research (09/2025)
Table 1: A summary of recovered PDEs from dataset using Adjoint and PDE-FIND method.
Problem (Ne, Na. Nagt..)
Method
Ex. Time [+]
Recovered PDE
Hest eq. (ID)
Adjoint
2.14
f = f.x + O(10-12)
(128,128)
PDE-FIND
2.66 0.02
ft far + 0(10-9)
Hest eq. (2D)
Adjoint
44.87 + 1.63
J: = Jaia: + Ja + ((10
(100,100,100)
PDE-FIND
796.01 + 11.91
f: Jz + far + -(10-9)
Burgers' eq. (LD)
Adjoint
3.37 0.24
J: = (12)2 + O
(128,128)
PDE-FIND
2.56 0.01
f 0.036f 0.094f 2ff. 0.002 0(10-1
Burgers' eq. (2D)
Adjoint
250.4 + 6.6
f: = (f') (f9 la31, 0010
(100,100,100)
PDE-FIND
914.93 +: 8.32
J: 1. 1.998 J=, + 1.998 +
KS eq. (1D)
Adjoint
9.34 + 1.2
f = - 0.5 0.5 )1 0(10-5)
(61,256)
PDE-FIND
1.12 0.11
ft 0.5f s. 0.5 1.972 -
Random Walk (1D)
Adjoint
1,258 + 0.38
f + 1.025 0.465
+0(0-))))
(50,100)
PDE-FIND
0.17 0.09
ft + 0.798 0.15-4.z
+0(10-9)
Reaction Diffusion eqs, (2D)
Adjoint
905.32 + 14.66
0.102 + 0.20 + 0.3u 0.3r
0.1av²
0.20A0 + 0.4a²
0110 ")
(200,70,70)
2 = 0.4% 0.33g + 0.2r 0.1.. -0.2we²
0.1. -0.6
+0,00
PDE-FIND
2231.54 31.52
M 0.10,
0. 0.200g + 0.3u 0.3r
0.1uv²
-020/' 0 0(10-))
R 0.4 + 0.30% 0.2v
0.1
0.2up
-.3. 0.1 0010
3.1.1
Heat equation
As a first example, let us consider measured data collected from the solution to the heat equation, i.c.
0f D 82
= 0,
(10)
Ōt
dr
with D D= 1. The data is constructed using the Finite Difference method with initial condition fxx 0) 5sin(2mx)<r(z - L) and mesh with N₂ 100 nodes in covering the domain Ω [0, L] with = 1 and N₂ 100 steps in with final time T = N, where At 0.05A |ID| is the step size and Ar L/Ng is the mesh size in x.
0.0
+ 0.-2,001
10
A(D]
10
-0.
||o|=")||
Coettikknt
10
10-7
Time (Ga))
10a
Ernor
10-1
104
-0.8
10-32
-10
100
20
20 Epochs
40
(100. 100)
(500, 100)
[1000 1000
11000 1000)
2001
(500, 1000
(1000 100)
11000 1000)
Epochs
(N,
(a)
(b)
(c)
(d)
Figure 1: The estimated coefficient corresponding to D in heat equation (a) and the -iomm error of all considered coefficients (b) using the proposed Adjoint method with N, = 100 and N, = 100. Also, we show Li-norm error of the estimated coefficients (c) and the execution time (d) using the proposed Adjoint method (blue) and FIN method (red), given data on a grid with N, E (100, 500, 1000) steps in t, and N, E (100, 1000) nodes in x.
We consider a system consisting of a single PDE (i.. =, = f, and p p) with one-dimensional input, i.c. n and x = z c IR, and d = d E N. In order to construct a general forward model, here we consider derivatives and polynomials with indices d, p e {1, 2, 3) as the initial guess for the forward model. This leads to 9 terms with unknown coefficients that we find using the proposed adjoint method (an illustrative
7

---

Published in Transactions on Machine Learning Research (09/2025)
While expect recover the coefficient that corresponds to D, we expect all the other coefficients (denoted by a²) ) to become negligible. That is what we indeed observe in Fig. I where the error of the coefficient for each term is plotted against the number of epochs.
derivation of the candidate terms be found in Appendix G.1
Table 2:
Comparison between Adjoint, PDE-FIND, WSINDy, and PDE-LEARN in recovering 1D Heat equation from noise-free data.
(z))
Method
Ex. Time [s]
Recovered PDE
TPR
(128,128)
Adjoint
2.14 = 0.01
Ji fax (110-12)
1
PDE-FIND
2.66 + 0.02
ft = faz + (
1
WSINDY
0.20 = 0.01
fr 10.777f.- 20.701
0.25
PDE LEARN
941.01 2
= + 0.002f + -0.009fee 0.006ff. f.007
0.22
-0.003(f=) 0.004 0.004f
Next, we compare the solution obtained via the adjoint method against PDE-FIND with STRidge optimization method. Here, we test both methods recovering the heat equation given data on the grid with discretization (Ne, Nz). As shown in Fig. the proposed adjoint method provides more accurate results across all data sizes.
We also point out that as the size of the data set increases, PDE-FIND with STRidge regression method becomes more expensive, e.g. one order of magnitude more expensive than the adjoint method for the data on a grid size (N,, = (1000, 1000). Table 2 we also compare the adjoint method with more recent methods such as WeakSINDy and PDE-LEARN, where, our surprise, PDE-FIND remains the strongest alternative that justifies its the baseline method here.
Next, consider the Heat equation in 2D, i.e.
0f 01
0
+ 0 =0 0
+
D
(11)
Ou7
Onl
with the initial condition
f(x,0) = exp(-b(pz x²)²) cos))oTT(2 (2: xx,)2)
(12)
where x. (0.5, 0.5)* and coefficients b = 20, inside the domain æ E [0. 112, as depicted in Fig. Again, have tested the adjoint method against other methods for variety of mesh sizes in Table Overall, the adjoint method seems to provide more efficient solution at larger data sets and higher dinensions.
200
1.
2.0
1.0
08-
4.8
04
0.8
0.5
08
2.6
06
0.6
0.2
6.4
02
12
0.2
-02
0.6
10
0.
0.4
[-
0.6
0.8
1.0
0.6
(a)
(b)
(c)
(d)
Figure 2: The initial condition described in Eq. 12 on 100 x 100 grid (a) and the evolution of the solution to the 2D heat equation after 20, 40, and 60 time steps (b)-(d).

---

Published in Transactions on Machine Learning Research (09/2025)
Table 3: Recovering 2D Heat equation using the Adjoint and PDE-FIND method given foese-free dataset for range of discrctizations. "The method failed in providing any solutions.
(No N,,,
Method
Ex.
Time [8]
Recovered PDE
TPR
(20,50,50)
Adjoint
11.4
Jt = SaYz J2,,, ((1055)
PDE-FIND
15.31
fur fegor 0(10-9)
WSINDy
120,100,100)
Adjoint
21.96
f1 fam + feim 0(10 5)
PDE-FIND
75.27
J1 = SELEI + +22m + 0(10-5)
WSINDY
0.139
ft -0.0009 fayaa + +0(109))
0.25
1100,100,100)
Adjoint
44.87
fi ferzi + s220 + O(10 6)
PDE-FIND
796.83
= faym feum 0(10-6)
WSINDy
0.152
ft -0.02166 f222222 + 0(1009)
0.25
3.1.2 Burgers' equation
As a nonlinear test case, let us consider the data from Burgers' equation given by
af
a(Af2)
(13)
Ōt
Ōr
where A -1. The data obtained with similar simulation setup as for heat equation (Section 3.1.1 for the time step, i.c. At 0.05Az/|A|
except
Similar to Section 3.1.1 adopt system of one PDE with one-dimensional input. We also consider derivatives and polynomials with indices d, E (1, in the construction of the forward model. This leads to 9 terms whose coefficients we find using the proposed adjoint method. As shown in Fig, 3 the proposed adjoint method finds the correct coefficients, i.e. 0d=11p== that corresponds to as well as all the irrelevant. ones denoted by a", up to machine accuracy in ((10) epochs.
0.0
0d-1pp2
10
-0.2
||eta Illls
Confiicient
10-2_
-0.4
0.6
10-2
10
-1.0
20
40
60
20
40
60
Epochs
Epochs
Figure 3: The estimated cocfficient corresponding to A (left) and the Lı-norm error of all considered coefficients (right) given the discretized data of 1D Burgers' equation during training.
Next, we compare the solution obtained from the adjoint method the one from PDE-FIND using STRidge optimization method. Here, we compare the error on coefficients and computational time between the adjoint and PDE-FIND by repeating the task for the data set with increasing size, i.c. (N,, N₂) e {(100, 100), (1000, 100), (1000, 1000)}. As depicted in Fig, 4 the adjoint method provides with more accurate solution across the different mesh sizes. Regarding the computational cost, while PDE-FIND seems faster on smaller data sets, as the size of the data grows, it becomes increasingly more expensive than the adjoint method. Similar to the heat equation, for the mesh size N2, = (1000, 1000) we obtain one order of magnitude speed-up compared PDE-FIND. Again, we compare the adjoint method WeakSINDy and PDE-LEARN as more recent alternative methods. As shown Table our surprise, PDE-FIND remains the strongest alternative that justifies its use the baseline method here.
9

---

Published in Transactions on Machine Learning Research (09/2023)
Adjoint
Adjoint
10
10a
10-3
Time.
10
101
109
[100, 100
(1000, 300
1000)
(100, 100)
(1000) 1000
(1000, 1000)
(1000)
(N,
Oh, (.)
Figure 4: Lı-norm enror of the estimated coefficients (left) and the execution time (right) for discovering the Burgers' equation equation using the Adjoint method (blue) and PDE-FIND method (red), given data grid with N, {100, 10003 steps in t, and N₂ E {100, 1000} nodes in x.
Table 4: Comparison between the Adjoint, PDE-FIND, WSINDy, and PDE-LEARN in recovering the ID Burgers' equation from noise-free data.
(N, N.)
Method
Ex. Time [s]
Recovered PDE
TPR
(128,128)
Adjoint
3.37 + 0.15
ft = (f²) 0(10-12)
PDE-FIND
2.56 + 0.02
2f 0.03 + f 0(10-4)
0.5
WSINDy
0.21 0.05
f. -00003..... + +00116ff
+ 0
0.25
PDE-LEARN
889.15 + 6.23
ft = 0.08f2
0.10 0.01 + +.03ff.++0066/ 0.034ffaz
0.2
+0.06f2fa- 0.05ff 0.04ffmmm
Next, consider the Burgers' equation in 2D, i.e.
0S Øt
0f2 05
+ A
+
(14)
Oz,
Ox
with and the initial condition
f(x. 0) = expl-b(a x²
(15)
with x. = (0.5, 0.5)T and coefficients 30, inside the domain æ E [0, 1/2, as depicted in Fig. Again, we have tested the adjoint method against the PDE-FIND method in variety of mesh sizes as shown in Table We observe that the adjoint method remains computationally advantageous at larger data sets.
t-
t-
0.
0.6
0.6
0.5
0.4
0.4
0.2
0.4
0.6
0.8
0.4
0.6
Figure 5: Initial condition (left) described in Eq. 15 and the solution of 2D Burgers' equation at final time (right).
10

---

Published in Transactions on Machine Learning Research (09/2025)
Table 5: Recovering 2D Burgers' equation using the Adjoint method without averaging Algorithm [ and PDE-FIND method given noiso-frec dataset for range of discretizations. "The WSINDy method failed in providing any solutions.
Method
Ex. Time [s]
Recovered PDE
TPR
520,50,50)
Adjoint
148.4
J: - 0.979( f2),, 0.988( ), 0(10-1)
PDE-FIND
15.23
f: 2.03f 2.02 fr. 0.027f fr.
0.029f
0.37
-0.035f 0.002fa, 0.002f,, 0(10-1)
WSINDY
(20,100,100)
Adjoint
145.6
Ft 0.989 )2 + 0.080( 0.004 (
0.42
+0.00M 0.107( ), -(10
PDE-FIND
85.34
ft -..01fff + 2.01 fa, 0.027 Sa, 0.012f2f.0
0.012 Sx,
0.3
-0.004 faa f3 0.002fz. - 0.002
0(10
WSINDy
0.1346
ft -0.001722 (f 1n2egu 0(10~9)
0.25
1000,100,100)
Adjoint
250.4
f = (f*),.. (f) 0(10
PDE-FIND
914.93
4 = 1.998 1.998
0(10-))
WSINDY
0.156
fe 00.1106271 ) + 0(10-
0.25
3.1.3
Kuramoto Sivashinsky equation
As more challenging test. case, let us consider the recovery of the Kuramoto-Sivashinsky (KS) equation given by
0f
05 Ōr
8f Or
8 0 Or
+ A
+B
+
(16)
ŌL
where A = -1, B 0.5 and = -0.5. The data is generated the grid (Ne, Nz) (64, 256) and the time step size At 0.01Ax*/11 /|C).
similar way previous sections except for
Here again, we adopt system of one PDE with one-dimensional input. As guess for the forward model, we consider terms consisting of derivatives with indices d E 2,, 4} and polynomials with indices b E {1, 21, leading to 8 terms whose coefficients we find using the proposed adjoint method. As shown Fig. the adjoint method finds the coefficient with error of ((10-5), yet achieving machine accuracy seems not possible.
Go 2
Ra-3,022
2.pe
m[A)
A(C)
Bp-2.p33
e(A)
eCC)
cLB)
|eee I11::
e(E)
|ele ")h
10"2
Coefficient
0.5
Error
10-5
Coeffidient
0.5
Error
10:
0.0
10-9-
0.0
10-0
10
10-11_
10-44
-1.0
-1.0
100
200
100 Epochs
200
100
200
100
200
Epochs
Epochs
Epochs
(a)
(b)
(c)
(d)
Figure 6: The estimated coefficients corresponding A, B, and the In-norm error all considered coefficients given the discretized data of the KS equation during training without (a-b) and with (c-d) active thresholding for Epochs Nuhr 100.
11

---

Published in Transactions on Machine Learning Research (09/2023)
Adjoint
10
Time (sec)
Enror
10-3
10-4
164,2500
(128,512)
186610244
(256.1024)
164,256)
(128.512)
(36,1024)
04. No)
(Ns, N₂)
Figure Li-norm eиror of the estimated coefficients (left) and the execution time (right) for discovering the KS equation using the Adjoint method (blue) and PDE-FIND method (red), given data on a grid with N₂ (64, 128, 256} steps in t, N₂ {256, 512, 1024) nodes in æ.
Table 6: Comparison between Adjoint, PIN-FIND, WSINDy, and PDE-LEARN method in recovering the Kuramoto Sivashinsky equation from data.
(N, N.)
Method
Ex. Time [s]
Recovered PDE
TPR
(64,256)
Adjoint
79.14 + 2.31
ft = -0.5f + 0.5f + (f 0(10~8)
1
PDE-FIND
26.20 +: 3.11
0.5fez + 0.5fmamx Fawr + 1.922 + 0.012 0(10-))
0.8
WSINDY
0.20 + 0.01
4 0.255 Jazz
0.2
PDE-LEARN
892.08 7.72
0.005 0.03 0.05fz + 0.02 0.05 ZMZ 0.04
0.26
00.07 + 0.03 0.04 J 0.05 faz + 0.08/ Jaza
-0.04 fazx f.07ffazzz 0.01 fazzz
Again, in Fig. 7 we make a comparison between the predicted PDE using the adjoint method against. PDE-FIND. particular, we consider a data set on temporal/spatial mesh of size (Ne, N₂) = (664, 256), (128,512), (256, 1024)) and compare how the error and computational cost. vary.
Similar to previous sections, the error is reported by comparing the obtained coefficients against the coefficients of the exact PDE in Lı-norm. Interestingly, the PDE-FIND method has 3 to 4 orders of magnitude larger error compared to the adjoint method. Also, in terms of cost, the training time for PDE-FIND seems to grow at a higher rate than the adjoint method as the (data) mesh size increases. Again, we made further comparison with more recent methods such as WeakSINDy and PDE-LEARN, As shown Table PDE-FIND remains the strongest alternative which justifies its the baseline method here.
3.1.4 Random Walk
Next, let us consider the recovery of the governing equation on probability density function (PDF) given samples of its underlying stochastic process. As an example, we consider the Itô process
dX = Adt + /2DdW
(17)
where A=1 1 is drift and = 0.5 is the diffusion coefficient, and W denotes the standard Wiener process with Var(dW)> = At. We generate the data set by simulating the random walk using Muler-Maruyama scheme
starting from X(t 0) for N₂ = 50 steps with time step size of At 0.01. We estimate the PDF using histogram with N₂ 100 bins and N. 1000 samples.
12

---

Published in Transactions on Machine Learning Research (09/2023)
Bs-1p=1
m(A)
||ela 1|
e(D]
1.0
1.00
0.8
10
0.4
10-4
0.2
0.0
100
200
100
200
Epochs
Epochs
Figure 8: The estimated coefficients corresponding to A and D (left) and the In-momm error of all considered coefficients (right) of the Bokker-Plancck equation as the governing law for the PDF associated with the random walk during training.
Adjoint
PDE-FIND
Adjoint
PDE-FIND
10"
10a
Time (sec)
109
Emrr
10-4
10-1
10-2
cat-oe (s)
Gat os)
Goo 1o0r'0or oor oot)
100, 109)
(50, 100.109
por-oor o01)
2000 001)
001)
(N.
Nal
Figure 9: Li-norm error of the estimated coefficients (left) and the execution time (right) for discovering the Fokker-Planck equation using the proposed Adjoint method (blue) and PDE-FIND method (red), given samples of its underlying stochastic process with N₂ E (50, 100} steps in t, N₂ = 100 histogram bins, and N₂ E (109, 104} samples.
Let us denote the distribution of X by f. Ttô's lemma gives the Fokker-Planck equation
0f
+ A 05
D 021 = Oz
(18)
Ōt Oz
Given the data set for f on a mesh of size (Ne, Na), can use Finite Difference compute the contributions from derivatives of in the governing law. Since this is one of the challenging test cases due noise, here we only consider three possible terms in the forward model, consisting of derivatives with indices d E (1,2,3) and polynomial power = 1. In Fig. S we show how the error of finding the correct coefficients evolves during training for the adjoint mot Clearly, the adjoint method seems to recover the true PDE with L error of O11--2) in its coefficients.
In Fig. g we make a comparison with PDE-FIND for the same number and order of terms as the initial guess for the PDE. We compare the two methods for range of grid and sample sizes, i.e. N, E (50, 100), N, 100, and N. E {103, 104). It turns out that the proposed adjoint method overall provides more accurate estimate of the coofficients than PDE-FIND, though at a higher cost. In Table 7 we show the discovered PDEs for both methods across the different iissretizations, and Table S shows further comparison wiih WSINDy and PDE-LEARN.
13

---

Published in Transactions on Machine Learning Research (09/2025)
Table Recovery of the Fokker-Plance equation, ie. ft+ 0.5 0, using the proposed adjoint method against PDE-FIND method given samples of the underlying stochastic process for various discretization
parameters.
N
2₂
N₂
Method
Recovered PDE
TPR
50
100
1000
Adjoint
St + 1.025/ - 0.465fax
PDE-FIND
fe + 0.798f. 0.454ffr
10000
Adjoint
fe 1.022f=
0.495ffx
PDE-FIND
La 0.818fx
0.496 0
1
100
100
1000
Adjoint
St + 1.010fz
0.543
0
PDE-FIND
ft + 0.863f
0.560
0
1
10000
Adjoint
ft + 1.015fz
0.589
0
1
PDE-FIND
Se + 0.894f. 0.612f
Table 8: Comparison between Adjoint, PDE-FIND, WSINDy, and PDE-LEARN in recovering the Fokker- Planck equation corresponding to the Random Walk from data.
(Ne,
Method
Ex. Time [s]
Recovered PDE
TPR
150,100)
Adjoint.
1.25 + 0.01
St + 1.025fz 0..655 fm 0)
1
PDE-FIND
0.17 + 0.02
ft + 0.798f= - -0.454f 0
1
WSINDy
0.19 +: 0.01
St = 7.87666
0.25
LEARN
69.10 = 2.13
= 0.0479ff + 0.0228f 0.0014f
0.75
3.1.5
Reaction Diffusion System of Equations
In order show scalability and accuracy of the adjoint method for a system of PDEs in higher dimensional space, let. consider system of PDFs given by
Ou
R" v) 0,
(19)
0t
2
dv
+
o
+
c
22
[b]
R² (", v) v)= =0
(20)
0i
where
(" v) = cau u + ch . + ouu2 + c5 uu'2v cu
(21)
(2, cyv + ciOu vuu + ch duu
(22)
We construct the data set by solving the system of PDEs Eqs. 19 20 using a 2nd order Finite Difference scheme with initial values
Up 4ww
3rZz
24) (L202 - 2 [L1
cos
Uo cos
sin
(L131 - xt) (L202 2²
La
where a = 100, and the coefficients
" 6 = 0001,,0022200033000000000,,00,,0
c" [ = 0.4, -0.3, -0.3,-0.2. 002,-0.1, 0.3, 0.2,
We generate data by solving the system of PDEs Eqs.
19
20
using the Finite Difference method and forward
Fuler scheme for N, 25 steps with time step size of At.
10-6 and in the domain Ω |0, [.] [0, T-2]
14

---

Published in Transactions on Machine Learning Research (09/2025)
u at t=
at t= T
1.0
5
5
2
2
X 0.5
0
0
-2
2
-5
-5
0.0
0.0
0.5
1.0
0.0
0.5
1.0
X1
X1
Figure 10: Solution the reaction diffusion system of PDEs time l T for a (left) and v (right).
e(cg) steg)
---
eicf) etey) etc])
e(cs)
---
e(e@)
ete()
steg)
etrp)
e(e!)
e(eg)
e(c))
10a
10-3
10-5
etc'l
10
10-9
10-2
10-13
10-12
100
200
100
200
Epochs
Epochs
Figure 11: Li-norm error in the estimated coefficients of the reaction diffusion system of PDEs during training.
where L L2 = 1 which is discretized using uniform grid with Na, N22 50 nodes leading to mesh size Ar1 Ar 0.02. In Fig. 10 show the solution to the system at time T N, for and v.
We consider a system consisting of two PDEs, i.e. N = dim(f) dim(p) = 2, with two- mio-dimensiona input, i.e. n = dim(x) = dim(d) = 2. Here, dfm(f) dim(Ima(ff)) where Ima(-) denotes the image (or output) of a function.
In order to use the developed adjoint method, we construct guess forward system of PDEs (or forward model) using derivatives up to 2nd order and polynomials of up to 3rd order. That is, dmax and Pmax 3. This leads to 90 terms whose coefficients we find using the proposed adjoint method (an illustrative derivation of the candidate terms can be found in Appendix G.2 The solution to the constructed model f R [u, u] as well as the adjoint equation for X is found using the samne discretization as the data set.
PDE on 4
10-7
||eaa Mis
10-5
10-7
10-9
11
100
200
Epochs
Figure 12: Li-norm error in the estimated coefficients of the irrelevant terms compared to the true reaction diffusion system of PDEs during training, i.e. [|e(("))] = ||a" ll--
15

---

Published in Transactions on Machine Learning Research (09/2025)
As shown in Fig. 11 the adjoint method finds the correct equations with error up (110-12) Furthermore. the coefficients corresponding to the irrelevant terms a tend to zero with error of O(10-11,, see Fig. 12 Furthermore, we have compared the adjoint method against PDE-FIND for a range of grid sizes in Fig. 13 We observe that the cost of PDE-FIND grows with higher rate than adjoint method as the size of the data set increases.
Adjaint
FDE-FIND
Adjoint
8000
10
6000
4000
2000
10
(25,50,50)
150,50.50)
(100,50.50)
1100,70,70)
(200.70.70)
(100,100,1000
100.10001000
1605000
100.0000
700.000000
100)
10000 100)
(200,100.00
12333
(100)
000700
100.
220
1ood
2002)
(N, N2)
(l. Nay Ne)
Figure 13:
Comparing the error and execution time the adjoint method (blue) to PDE-FIND method (red)
against the size of the data set for the tolerance of 10-7 in the discovered coefficients.
3.1.6 Wave equation
Consider wave equation
f(x,t) sin(x - t)
(23)
which is a solution infinite PDEs. For example, one class of PDEs with solution f
+ + (I: 1)5zza + c(Jan Sarna) Vk and c e R,
(24)
defined a domain æ E [0, 2 and T = 1. We create a data set using the exact f on a grid with N; = time intervals and N₂ 100 spatial discrctization points.
Let us consider similar setup the heat equation example als indices d {1,..., 6} and p = 1 as the initial guess for the forward model. This leads to 6 terms with unknown coefficients a. Here, we enable averaging and use a finer discretization in time (100 steps for forward and backward solvers in each time interval) to cope with the instabilities of the Finite Difference solver due to the inclusion of the high-order derivatives. We also disable thresholding except at the end of the algorithm.
with derivatives and polynomi-
As shown in Fig. 14 the proposed Adjoint method returns the solution
+ 0.996f.
(25)
which is the PDE with the least number of terms compared all possible PDEs. We note that for the same problem selting, PDE-FIND identifies the same form of the PDE, i.e.
f 0.9897 0
(26)
16

---

Published in Transactions on Machine Learning Research (09/2025)
1.0
0.5
Coefficient
0.8
Rg-2p=1
R/-1p=1
0.6
0jeipml
0.4
Rjaapmi
-0.5
0.2
0csSpe1
0ga0p=1
-1.0
2
200
400
Epochs
Figure 14: Profile of f at = 0 and t = 1 (left) and the evolution of considered coefficients during adjoint optimization (right)
3.2 Partial observations in time
Here, we investigate how the error of the discovery task increases when only a subset the fine data set is available. Consider the heat equation presented in section 3.1.1 and consider a data sot created by solving the exact PDE using the Finite Difference method with A T/N, where N, = 1000 and Ar L/Ng and N₂ 1000.
Let us assume that we are only provided with a subset of this data sct. As test, let us take every v time step as the input for the PDE discovery lask, where v E (1, 4, 8, 16]. This corresponds to using {100, 50, 25, 12.5, 6.25]% the total data set. By doing so, the accuracy of the Finite Difference method in estimating the time derivatives using the available data deteriorates, leading to large error PDE discover task for PDE-FIND method.
However, the adjoint method can use a finer mesh in time compared to the data set in computing the forward and backward equations and only compare the solution to the data on the coarse mesh where data is available. We use N, = 1000 for the forward and backward solvers in the adjoint method, and impose the final time condition where data is available. As shown in Table g and Fig. 15 the proposed adjoint method able to recover the exact PDE regardless of how sparse the data set in time.
We emphasize that while adjoint method can use finer discretization in time than the one for data G in solving forward and backward equations, it is bound to similar coarser spatial discretization G. This is due to the fact that the data points f are used for the initial condition of the forward model eq. g and the final condition of the backward adjoint equation eq.
Adjint
Adont
10-2
e(D1
e(D)
205
|ete") b
eeee 1:
10-3
10
10-2
Time (ooc)
10a
10-3
10-9
10-4
10-4
10-42
9.
10a
40
25
50
75
100%
23%
12.5%
6,25%
100%
50%
23
17.5%
0.29 %
Epochs
Epochs
(a)
(b)
(c)
(d)
Figure 15: Evolution the Li-norm error coefficients considered terms using adjoint method when only (a) 50% and (b) 6.25% (b) of the data available. Error and execution time of the adjoint method (bluc) and PDE-FIND method (red) in finding the coefficients of true heat equation given sparse data set in time in (c) and (d).
17

---

Published in Transactions on Machine Learning Research (09/2025)
Table 9: Recovering the heat equation given sparse dataset in time. Here we rounded the coefficients up three decimals.
Meshod
Recovered PDE
TPR
100
Adjoint
f2
PDE-FIND
faz
50
Adjoint
fax -0
PDE-FIND
f. 0.900 + 0.177f - 0.261 0.089ff2
0.25
-0.011 [. 0.0033 fz= 0.001f
Adjoint
ft faz
PIN-
0.532f 0.778 f 0.268ff. 0.010 fax 0.035f f,-
0.999
0.25
0.003 fuee
12.5
Adjoint
fax
PDE-FIND
f. 0.908f
1.264f
f.803f 0.638
0
0.22
0.007 f feo- -0.026 f.x
0.001 f..=
Adjoint
faz
FIND
J. 0.999/ 2.769
f.051
1.398 11,
0.185 3
0.22
-0.055/55..
0.016 J...
0.002 Jers
FIND
FIND
avg.
aay.
102
10*
Time (soc)
Time (sec)
10
Eraor
200
10-4
10
10a
0.001
0.093
0.11
0,001
0.001
0.005
0.01
0.01
Naise 195)
Noise (16)
Noise (%)
Noise (%)
(a)
(b)
(c)
(d)
Figure 16: Error and execution time of the adjoint method with (green) and without averaging the gradients (blue), along with the PDE-FIND method (red) in finding the coefficients of the true PDE, i.e. heat equation (a)-(c) and Burgers' equation (b)-(d), given noisy data.
3.3 Sensitivity to noise
Let us investigate how the error increases once noise is added to the data set. In particular, we add noise to each point of the data set for f" via f* c), where the noise to signal ratio is (cf*)/f" c ~ NNO0 o2) with o²) denoting a normal distribution with zero mean and standard deviation of a. As test cases, revisit the heat (section 3.1.1] and Burgers' equations (section 3.1.2 with added noise of with e {0.001, 0.005, 0.01, 0.1} %. Before searching for the PDE, first denoise the data set using Singular Value Decomposition and drop out terms with singular value below a threshold of O(10-4). Here, we report the expected coefficients of discovered PDE and execution time by repeating the process 10 times.
As shown Fig. 16 adding noise to the dataset deteriorates the accuracy finding the correct coefficients of the underlying PDE for both the adjoint and PDE method. We observe that the adjoint method, both with and without gradient averaging, less susceptible to noise compared to PDE-FIND, albeit at higher computational cost. Additionally, averaging the gradients in the adjoint method improves the accuracy around two orders of magnitude at higher computational cost.
Let. us again revisit the Heat and Burgers' equation test. cases in 2D, and compare the adjoint method with averaging gradient algorithm 2 against PDE-FIND method without applying any noise reduction Table 10- II Although the accuracy of both methods deteriorates with noise, the adjoint method seems provide more reliable solution.
18

---

Published in Transactions on Machine Learning Research (09/2025)
Table 10: Recovering 2D Heat equation given noisy data on the grid N,,, = (50, 100. 100) using the adjoint method with averaging and learning rate parameter / = 0.02 in Algorithm 2 and PDF-FIND. Given the raw noisy data, PDE-FIND was not able to find any PDE, as it lead to run-time errors. Here, we also report the discovered PDE for denoised data with SVD.
0
Denoised
Method
Time [s]
Recovered PDE
TPR
0.1 %
Adjoint
108.5
ft = fr,,, fann + 0.002 0.003
f: SamZ + SNAZ, + 0(10 2302 + O(10-t)
0.6
yes
110.5
PDE-FIND
235.40
Si 1.003 2m + 1.02 2o0e 14.44 37.70 + f2 f0.14 [am f f2 Jagus 0.15 f.95 0,10-
0.33
yes
227.1
S = Jan + fears + 0(10-9)
1 %
Adjoint
110.2
fr 0.963 + 1.008f 0.147
0.5
+0.010f + 0.002 + 0(10
yes 109.19
ft = faje, + fazaz + 0.022f. + 0.027 + 0.031f
2222 +
+1100-9)
0.5
PIE-FIND
yos 245.3
ft 0.36 + 1.01fa + 1.01fag + 57.08 + 110.54f 32.68
0.18
+1.78f f., + 3.20 f.z 1.06f 0.77f + 0.20 faym
+0.44f fay 0.100 o, + f.2 + 0.64ff 0(10-2)
5 %
no
Adjoint
105.43
ft 1.30 11 + 1.322 - 0.135 + 0.50 + +.18f - -0.23f2
0.3
-0.19/2 - 0.22 0.20 + 0
yes 106.21
fe = 1.06
11,07 0.91 + + 533, + +0.12 0.11.2.
0.3
-0.15f2, + f..11
0.10 +0(10-
PDE-FIND
yos
248.53
4- 1771.33/ 2369.4992 -971133/-2300449f2 $ - 11.41
0.16
-20 -66446 4,40/f2, + T,77f 147 + 7 7.40
-0.81 +034f + 1.45 + 1.46 -199 fro
+7.37 + +281 0(10-))
Table 11: Recovering 2D Burgers' equation given raw noisy data on the grid (Ne, No,, No,) = (50, 100, 100) using the adjoint method with averaging and learning rale parameter B = 0.01 Algorithm and PDE-FIND. Given the raw noisy data, PDE-FIND was not able to find any PDF, as it lead to run-time errors. Here, we also report the discovered PDE given noisy dataset that is denoised using SVD.
Denoised
Method
Time [s]
Recovered PDE
TPR
0.1 %
Adjoint
210.2
f. 0.971(f)) + ( )x + 0.024f2, + 0.027 + 0,10 9)
0.6
yes
214.09
== + 0.992( +0,, 0(10-3)
1
PIE-FIND
280.71
0.115
42.29/ + 140.06f 65.644 + 0.21fff + f.51 - f.
0.17
+0.17ffa, + #4.22f - 2.14
0.77
0.94ff faLa, -
0.31 fum
+0.79/
1.01/ / + 0.38
yes
99.32
fe == 1.900.f.. 1.990
+0(00-
1 %
Adjoint
220.2
fi f ), + 0.511
) 0.09 0.51())) 0.000""__000000000000000......."
0.09f
0.27
+0.13 (.01(f") 0.13
0.00( 0(1003)
yes
231.0
f= 0=965))) 0...8. +.00, +0.018ff. 0.us2(f),,, -0.0181
0.12
PDE FIND
279.41
=== 0.18 0.18 +.19 -
-3.1 Ja
0.06
0.75f 17.51 31 +0
+3.21/ fz, 0.78/
+ -
1128f
yes
290.01
+ 32.98 Je = 14.53
- 1511f - f
0.14
+0.70f + 0.66 fa.,, 0.69f
O(10-)) t
3 %
Adjoint
240.9
f- 0.51 0.315 + ",,0 0.788 -,7
033
-0.788 | -0000
yo8
242.90
f- - 0,300, 0500 0.155/)))_,_,,111,,, 0004,,
0.6
no
PDE-PIND
yes
288.91
ft = -= +324 1.23/
0.07
-1.22 0.19 0.18 0.78ff 0 f
+0.75/fo0m 0..6f 0(10
19

---

Published in Transactions on Machine Learning Research (09/2025)
3.4 Addressing ill- fllpposedness
There may exist more than one PDE that replicates the data set. Therefore, the PDE discovery lask is ill-posed due to the lack of uniqueness in the solution. This is an indication that further physically motivated constraints are nceded to narrow the search space to find the desired PDE. However, among all possible PDEs, which PDE is found by the Adjoint method with the loss function defined as Eq. 3
To answer this question, let us consider a simple example of the wave equation
f(x,t)= = sin(x -t)
(27)
which is a solution infinite PDEs. For example, one class of PDEs with solution f ft + kfz + (k- - 1)fezz + clfex + fezzz)= = E CER
(28)
defined in a domain x E [0, 2n] and T 1. We create a data set using the exact f on grid with N, 10 time intervals and N₂ 100 spatial discretization points. Let us consider similar setup the heat equation example 3.1.1 with derivatives and polynomials indices d E f1.... 6} and p p= as the initial guess for the forward model. This leads to 6 terms with unknown coefficients
Here, we enable averaging and use a finer discretization in time (100 steps for forward and backward solvers in each time interval) to cope with the instabilitics of the Finite Difference solver due to the inclusion of the high-order derivatives. We also disable thresholding except at the end of the algorithm.
The proposed Adjoint method returns the solution
ft + 0.996f = 0
(29)
which is the PDE with the least number of terms compared to all possible PDEs. We note that for the same problem setting, PDE-FIND finds
Je - 0.9897f = =
(30)
The identified form of PDE be explained by the use of regularization term in the cost function 3 which enforces the minimization of the PDE coefficients. Clearly, the regularization term may be changed to find other possible solutions of this ill-posed problem.
3.5 Incomplete guessed PDE space
In this section, we investigate the outcome of the adjoint method when the exact terms not included in the initial guessed PDE form. Here, we define the space of PDE where the exact terms are included in the general forward model [ as complete. If the considered general form of PDE I does not include all the terms of the exact PDE, we denote that as an incomplete guessed PDE space.
Let us take the data from the numerical solution to Burgers' equation used in section 3.1.2 with discretization N; 22 100. For the complete forward model, we again consider derivatives and polynomials with indices d.p E (1, 3) in the construction of the forward model. This leads to 9 terms whose coefficients find using the proposed adjoint method. For the incomplete space of PDE, we take derivatives and polynomials with indices as d E (1,2,3} and E {1,3), leading to 6 terms. Clearly, the incomplete guessed PDE space does not include the term 0, 20 Now, we would like to see which PDE is returned by the adjoint method.
In Fig. we made a comparison between the evolution of coefficients and L2 norm error of the estimated forward model against the data. While the complete space monotonically converges to the exact solution up to machine accuracy, the incomplete space of PDE delivers another PDE, i.e.
Øf 8² +
Øt a. (0.04f 0.01f*) = 0,
(31)
with the relative L2 error of O(10-) between forward model estimation and the data points. The fact that the Ta error between and f does not decrease is an indication that the considered space of PDE is
20

---

Published in Transactions on Machine Learning Research (09/2025)
10-3
Ro-10=7
10-1
00-1p=2
""11:1 11s
10-1
00.2.001
0..
10
Gro2.-2
0.6
10:44
0g-2pmi
-0.8
10-17
10-16
20
Epochs
Epochs
Epochs
10"1
0.04
2.47 10
10-1
0.03
IIf "'llsl lls
00-3,03
2,465 10
10-1
0.02
0g-2.p=1
16-4
0c-2,pm3
0.01
02-2,po2
2.46x10
10-4
Bootpsl
2.455 10
10
1
-0.01
20
Epocts
Epochs
Epochs
Figure 17: The adjoint method applied to the Burgers' equation for complete (Lop) and incomplete (boltom) space of guessed PDEs.
incomplete and additional terms must be included. We note that here assumed there is noise in the data set. However, in the presence of noise, the L2 error between f and f may stagnate at the noise level, which makes the analysis on the completeness of the PDE space more challenging.
4 Discussion
Below we highlight and discuss strengths and weaknesses of the proposed adjoint method.
Strengths.
The proposed method has several strengths:
1. The proposed adjoint-based method of discovering PDEs can provides coefficients of the true governing equation with significant accuracy.
Since the gradient of the cost function with respect to parameters are derived analytically, the optimization problem converges fast. In particular, the adjoint method becomes cheaper than PDE-FIND as the size of the data set increases. The adjoint method by construction finds the optimal relation between the gradient of the cost function and the error in the data points. This was achieved by finding the extremum of the objective functional using the variational derivative.
We note that a clear difference from the point-wise loss - f used in the PDE-FIND method is that for the adjoint method weights the error at discrete points with the Lagrange multipliers; see Eq. g and the final condition Eq. 7 For Adjoint method more data often means finer resolutions which allows solving the forward/ prrward/backward backward equations more accurately on finer grid, and obtain higher accuracy in estimating the gradients that lead to less number of epochs in training, see Appendix
Since the adjoint method uses PDE solver find the underlying governing equation, there is a guarantee that the recovered PDE has solution and can be solved numerically.
The adjoint mothod can use a finer mesh in time compared to the available discretization of the data set. This allows an accurate recovery of the underlying PDE compared to the PDE-FIND, where the error the latter increases as the data set gets coarser since it estimates derivatives directly (either with Finite Difference a polynomial fit) using the given data set.
21

---

Published in Transactions on Machine Learning Research (09/2025)
Weaknesses. Our proposed method has some limitations:
In order to use the proposed adjoint method for discovering PDEs, a general solver of PDEs needs to be implemented. Here, used Finite Differences which can be replaced with more advance solvers. Clearly, the proposed adjoint method is most effective when there is prior knowledge of the underlying PDE form, and an appropriate numerical solver is deployed. We note that an inherent. limitation of the proposed adjoint method is the possibility of encountering either ill-posed forward or backward equations during optimization, which limits the time step size.
2. In this work, we used the same spatial discretization as the input data. If the spatial grid of input data is too coarse for the PDE solver, one has to use interpolation to estimate the data on a finer spatial grid that is more appropriate for the PDE solver.
3. In this work, we made the assumption that the underlying PDE can be solved numerically. This be limitation when there are no stable numerical methods to solve the true PDE. In this scenario, the proposed method may find another PDE that is solvable and fits to the data with a notable error.
4. Similar to PDE-FIND and similar symbolic regression methods, the Adjoint method considers a library of symbolic terms for the PDE. This can be a limitation when no prior information about the underlying dynamics is available.
5
Conclusion
In this work, we introduce novel mathematical method for the discovery of partial differential equations given data using the adjoint method. By formulating the optimization problem in the variational form using the method of Lagrange multipliers, we find an analytic expression for the gradient of the cost function with respect to the parameters of the PDE as a function of the Lagrange multipliers and the forward model estimate. Then, using variational calculus, we find a backwards-ing tekkwards-in-tinne time evolution equation for the Lagrange multipliers which incorporates the error with a source term (the adjoint equation). Hence, we can use the same solver for both forward model and the backward Lagrange equations. Here we used Finite Differences to estimate the spatial derivatives and forward Fuler for the time derivatives, which indeed can be replaced with more stable and advanced solvers.
We compared the proposed adjoint method against PDE-FIND in several test cases. While PDE-FIND seemss to be faster for small size problems, we observe that the adjoint method equipped with forward /bacwaard solvers becomes faster than PDE-FIND as the size of the data set increases. Also, the adjoint method can provide machine-accuracy in identifying and finding the coefficients when the data set is noise-free. Furthermore, in the case of discovering PDEs for PDFs given its samples, both methods seem suffer enormously form noiso/bian associated with the finite number of samples and the Finite Difference on histogram. This motivales the use of smooth and least biased density estto in these methods such Tohme et al. (2023) in future work. In the future work, we intend to combine the adjoint-based method for the discovery of PDE with PINNs as the solver instead of Finite Difference method. This would allow us to handle noisy and sparse data as well as deploying larger time steps in estimating the forward and backward solvers.
Acknowledgements
T. Tohme was supported by the MathWorks Engineering Fellowship. M. Sadr thanks Dr. Andreas Adelmann for his support and Prof. Thomas Antonsen for his constructive discussion on the adjoint method. This research was also supported by the Center for Complex Engineering Systems at the Massachusetts Institute of Technology (MIT) and the King Abdulaziz City for Science and Technology (KACST).
References
E Paulo Alves and Frederico Fiuza. Data-driven discovery of reduced plasma physics models from fully kinetic simulations. Physical Review Research, 4(3):033192, 2022.
22

---

Published in Transactions on Machine Learning Research (09/2025)
Thomas Antonsen, Elizabeth J Paul, and Matt Landreman. Adjoint approach to calculating shape gradients for three-dimennional magnetic confinement equilibria. Journal of Plasma Physics, 85(2):905550207, 2019.
Thomas Blumensath and Mike E Davies. Iterative hard thresholding for compressed sensing. Applied and computational harmonic analysis, 27(3))265-274 2009.
Josh Bongard and Hod Lipson. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104((4):9943 9948, 2007.
Steven Brunton, Joshua Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academyy of sciences, 113 (15):3932 3937, 2016.
Russel Caflisch, Denis Silantyev, and Yunan Yang. Adjoint dsmc for nonlinear boltzmann equation constrained optimization. Journal of Compulational Physics, 439:110404, 2021.
Kathleen Champion, Bethany Lusch, J Nathan Kutz, and Steven L Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 2019.
Zhao Chen, Yang Liu, and Hao Sun. Physics-informed learning of governing equations from scarce data. Nature communications, 11(16, 2021.
Bryan C Daniels and Ilya Nemenman. Automated adaptive inference of phenomenological dynamical models. Nature communications,, 6(1):8133, 2015a.
Bryan C Daniels and Ilya Nemenman. Efficient inference of parsimonious plocoeeooooaa models of cellular dynamics using s-systems and alternating regression. PloS one, 10(3):o0119821, 2015b.
Junming Duan and Jan S Hesthaven. Non-intrusive data-driven reduced-order modeling for time-dependent parametrized problems. Journal of Computational Physics, 497:112621, 2024.
Solve Fidnes and Kjetil Olsen Lye. Pseudo-hamiiltniian neural networks for learning partial differential equations. Journal of Compulational Physics, page 112738, 2024.
H Pearl Flath, Lucas C Wilcox, Volkan Akçelik, Judith Hill, Bart van Bloemen Waanders, and Omar Ghattas. Fast algorithins for bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial hessian approximations. SIAM Journal on Scientific Computing, 33(1):407-432, 2011.
Alcssandro Geraldini, Matt Landreman, and Elizabeth Paul. An adjoint method for determining the sensitivity of island size to magnetic field variations. Journal of Plasma Physics, 87(3):9058700) 2021.
Dimitrios Giannakis and Andrew J Majda. Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability. Proceedings of the National Academy of Sciences, 109(7): 2222 2227, 2012.
Mariano Giaquinta and Stefan Hildebrandt. Calculus of variations I, volume 310. Springer Berlin, Heidelberg, 2004.
Raul Gonzalez-Gareía, Ramiro Rico-Marthnez, and Toannis G Kevrekidis. Identification of distributed parameter systems: A neural net based approach. Computers & chemical engineering, 22:5965-8968, 1998.
Thomas JR Hughes, Gonzalo R Feijóo, Luca Mazzei, and Jean-Baptiste Quincy. The variational multiscale method paradigm for computational mechanics. Computer methods in applied mechanics and engineering, 166(1-2):3- 1998.
Antony Jameson. Acrodynamic shape optimization using the adjoint method. Lectures at the Von Karman Institute, Brussels, 2003.
Kadierdan Kaheman, J Nathan Kutz, and Steven Brunton. Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A, 476(22222):02200277) 2020.
23

---

Published in Transactions on Machine Learning Research (09/2025)
Ioannis G Kevrekidis, C William Gear, James M Hyman, Panagiotis G Kevrekidis, Olof Runborg, Constantinos Theodoropoulos, et al. Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis. Commun. Math. Sci, 1(4):715 162,7,1-7666--- 2003.
John H Lagergren, John T Nardini, G Michael Lavigne, Erica M Rutter, and Kevin B Flores. Learning partial differential equations for biological transport models from noisy spatio-temporal data. Proceedings 2020.
of the Royal Society A,
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In International conference machine learning, pages 3208 3216. PMLR, 2018.
Zichao Long, Yiping Lu, and Bin Dong, Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid deep network. Journal of Computational Physics, 399:108925, 2019.
Andrew J Majda, Christian Franzke, and Daan Crommelin. Normal forms for reduced stochastic climate models. Proceedings of the National Academy of Sciences, 106(10):3649 3653, 2009.
Niall M Mangan, Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Inferring biological networks by sparse identification of nonlinear dynamics. TEEE Transactions on Molecular, Biological and Mulli- Scale. Communications, 2(1):52 63, 2016.
Takashi Matsubara, Ai Tshikawa, and Takaharu Yaguchi. Deep energy-based modeling of discrete-time physics. Advances in Neural Information Processing Systems, 33:13100 13111, 2020.
Daniel A Messenger and David M Bortz. Weak sindy for partial differential equations. Journal of Computa- tional Physics, 443:110525, 2021.
Tgor Mezić. Analysis of fluid flows via spectral propertios of the koopman operator. Annual review of fluid mechanics, 45:357-378, 2013.
Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Compulational Physics, 357:125 141, 2018.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: A deep
learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-7000 707, 2019.
Anthony John Roberts. Model emergent dynamics in compler systems, volume 20. SIAM, 2014.
Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Data-driven discovery of partial differential equations. Science advances, 3(4) 1626144 2017.
Nihar Sawant, Boris Kramer, and Benjamin Peherstorfer. Physics-informed regularization and structure preservation for learning stable reduced models from data with operator inference. Computer Methods in Applied Mechanics and Engineering, 404:115836, 2023.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324 (5923):81 :8-22) 15, 2009.
Michael D Schmidt, Ravishankar R Vallabhajosyula, Jerry W Jenkins, Jonathan E. Hood, Abhishek S Soni, John P Wikswo, and Hod Lipson. Automated refinement and inference of analytical models for metabolic networks. Physical biology, 855)0055011, 2011.
Robert Stephany and Christopher Earls. Pde-learn: Using deep learning to discover partial differential equations from noisy, limited data. Neural Networks, 174:106242, 2024.
George Sugihara, Robert May, Hao Ye, Chih-hao Hsich, Ethan Deyle, Michael Fogarty, and Stephan Munch. Detecting causality in complex ecosystems. science, 338(6106):446 5006 2012.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Stolistical Methodology, 58(1):267-285 1996.
24

---

Published in Transactions on Machine Learning Research (09/2025)
Tony Tohme, Dehong Liu, and Kamal Youcef-Toumi. GSR: A generalized symbolic regression approach. Trnsaccions on Machine Learning Research, 2022.
Tony Tohme, Mohsen Sadr, Kamal Youcef- Toumi, and Nicolas G Hadjiconstanntion. Messy estimation: Maximum-entropy based stochastic and symbolic density estimation. arXiv preprint arXiv:2306. 04 120, 2023.
Henning Voss, Paul Kolodner, Markus Abel, and Jürgen Kurths. Amplitude equations from spatiotemporal binary-fluid convection data. Physical review letters, 83(1))3442,, 1999.
Hao Ye, Richard J Beamish, Sarah M Glaser, Sue CH Grant, Chih-hao Hsieh, Laura Richards, Jon T Schnule, and George Sugihara. Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling. Proceedings of the National Academy of Sciences, 112((3):E1569 E1576, 2015.
25

---

Published in Transactions on Machine Learning Research (09/2025)
A Derivation of the Adjoint equation
In this section, we provide a detailed derivation of the adjoint equations presented in this paper. Let us consider the forward model
Oef: + aidd V V =0 d.p
(32)
for i = I,....N. N. multipliers to formulate a cost functional for the function f that has the minimum error from the data points f" with the constraint that also solves the forward model. For simplicity, let us consider the cost functional only within each time interval [r60) x(+11), i.e.
order to find the parameters Oi,dpp of the such model, we use method of Lagrange
(li (li41l (((tt)) f((x11), (20+1))2++ / (x(x,,)C,,) ()ddd + col col|o
cffl Σ
(33)
Assuming the solution f and the Lagrange multipliers A are sufficiently smooth, we derive functional derivatives, more precisely Gateaux derivatives Giaquinta and Hildebrandi (2004), of C with respect f within the time interval [t(), (0+1).. We perform these operations first for the term denoted by
81[f] = lim
- [ff + edfieil
d
(34)
de
lim
d
Σ (f) (2(00, (042)) - fx( (()t)) - cff((z (6)+1))))
(35)
de
lim Σ (z) ( (f (,)0)),)++))))) - (z) ((0),,))))))) - c6f((z(1)) (000(())66-1))))) (((--1))))
(36)
(k) (f+ (x) ((+t1)) - (((x(())((-1111
(37)
and the second term denoted by J,
8J[ = 1im d
J[f + cōfie:))
(38)
c-0 de
d
/ > + +edsieildadd
lim c-0 de
(39)
-lim( d <0> de / 'x.011it + c65l + +> Σ
ai.d v(d) f8 ...(Ji p...)t.... (Kddd)))
(40)
d.p
= =lim( d <0 de [ -0l1-U +65 + Σ 1-1) |4| 04 dp
V [AA 1 f3- + Pt. S. dxdt
d.p
(0+1)
+ (Ai(fi + +effi))
(41)
de
= -limg /
-0611115f +
(-1)d ds p ... + 1 PPN (d) [> dxdt)
...
d.p
+ +AroSi) 16ty
(42)
(
= =/0f -0,1a + (-
fr ((d) [idl dxdt + (2,8f))
(541)
aidd
(43)
fi
165)
d.p
26

---

Published in Transactions on Machine Learning Research (09/2025)
= =f 85 -..(A.] + (-1))daa d.p & L/P]v/ P.] dadd + de(ar, (0+1))181166 80 10-1))
(44)
d.p
Here, e; denotes the standard basis vector in the direction of ith-dimension. In this derivation, we used integration by parts, the divergence theorem, and considered compact support for X in the spatial solution domain Ω i.e. + 0 on the boundaries An. Since we assume that the solution at the imitial time is given, there no variation of C at the beginning the time interval, that is, off(æ, (()) = Hence, the total variation of the cost functional with respect to the function f for the time interval (t19) e(++1) becomes
80|f1= = Σ
Σ
2(F((al0), - fi(xllh), (i
i=1
/
Odi
((d))(aadd (dl)of:dadt
òt
d.p
+ Σ ,
(45)
Now that have found the total variation of the cost functional C with respect to f. we can find the adjoint equation by considering only variation with respect to f in E and ignoring , leading
r0+1),
sCifl
=0
(46)
f(a,00+11)
/
0d:
+ Σ
1)ld
dj.d.p
8f,dadt 0
(47)
Öt
d.p
Σ d.p
1)
ag V [ [1,]
(48)
Ot
This equation is called adjoint equation presented in oquation 6 Similarly, we can find the final condition for the adjoint equation by only considering variation of f at final time = +0++1 and ignoring of, which we used to denote the variation of for E (t6), ((+11),, leading to
6c[f]
= (k) ((+1)) 2(ff (x) ((+t)) (()t)))))))))))) f)( (((+1))))
(49)
((x11) viE(US),
for = I,.., N and j = O,..., N₂
B Derivation of the Adjoint gradient with respect to PDE parameters Consider the cost functional
N c(fl = Σ (e)((( (66t)) (t) f)()())) / (x((x,)) (x) + colla
(50)
where
C.ff = Ō.fe + Σ ai,d Void [
(51)
27

---

Published in Transactions on Machine Learning Research (09/2025)
Assume that have solved the adjoint equation equation 6 and found the Lagrange multipliers À for each time interval (10), (0-11, Next, we can find the gradient of C with respect to PDE parameters simply by
OC
=
As [P + dd.
(52)
das.d.p
Using inlegration by parts, Divergence theorem, and the fact that À has a compact support on the boundaries, obtain
aC = V(di (TsPP (-1) f fr (PVid (2i]dadtt + 20.0rid.pp
(53)
da;
-(-1)ld) / fP V DA,ldædt + 2cn0riddpp
(54)
C Error in the numerical estimate of the Adjoint gradient
Consider the second-order central finite difference scheme as spatial and the first-order Euler as the temporal discretization scheme for the forward equation [ and backward equations equation Let us denote the discretized approximation of the solution with f and 2, leading
f = f + +00))++000)
(55)
\= - 0
(56)
where h is the spatial spacing and At is the time step size. By plugging the discretization into equation and using the same second-order numerical integration scheme in and the first-order scheme in it can be shown that
aC da;,d.p = (-1)ld) / ( (f++0())+ (v)) [M + O(h²) d(tt) + 20g0r,d.p
(57)
= (-1)ld) / f V + 2egas.d.p
(58)
OC
+0(0))+
(59)
doi.d.p
OC denotes the discretized estimate of the gradient of cost function. Therefore, the adjoint estimate where dog.d.p
used in this work sccond order in æ and first order in Clearly, this may be improved using stable higher-order scheme. Interestingly, the accuracy of the adjoint gradient is only function of mesh spacing and time step size, and not number of data points.
D Justification for the choice of the learning rate
In the proposed adjoint method, considered the update rule
ac
Qi,d + ai.d.p Dog.d.
(60)
for = 1,... . N. Here, we give a justification for our choice of the learning parameter From the expression for the gradient of cost function with respect to parameters 4 i.c.
28

---

Published in Transactions on Machine Learning Research (09/2025)
ac
--1)1a)) f pd PPV (i]daad + d.p
(61)
doi.dd.
we can see that
0C - - d0
dx) 0(V(
(62)
d.s. S 0h Id
(63)
where h min(Ar). So, the magnitude of the gradient scales exponentially with the order of the derivative d. The highest order terms, ie. the terms with d = dmax max(ld)), have the largest magnitude for their gradients. This means that by taking constant learning rate 7) the adjoint method would find the coefficients of the highest order terms first. This effect leads to the non-uniform convergence of the adjoint method.
In order to enforce uniform convergence on all PDE parameters, in this paper we consider
3 min(Ar) (Ar) dmax
(64)
as the learning rate which encodes the scaling with respect to the order of derivative for each PDE term. With this choice of learning rate, we have
OC dos
n|
5 - (mmx))) ((h 11111)
(65)
S ( damax +1 d1mx
(66)
for all i, d.p. Hence, our choice of n, i.e. Eq. 64 enforces uniform convergence on all PDE parameters.
E Further Details on Adjoint method
Here, we present flowcharts to illustrate the proposed adjoint algorithms 12 with and without averaging the gradients in Fig. 18 We also provide instructions to help practitioners to deploy the Adjoint method for their data set.
Diverging Coefficients. There are two scenarios in which the learning of coefficients in the Adjoint method can become divergent. First, the intermediate guessed PDE may become so unstable that numerical solution to either forward or backward adjoint equations diverges even for short time intervals E [(0+)), {6)]. Second, the time step size from the dataset {6++) {0) is too large for the underlying numerical solver, leading to diverging coefficients. coefficients during training, one should reduce the time step size for the forward and backward Adjoint equations which are independent of time step size in the dataset. Second, in these scenarios, we recommend deploying smaller learning rate 1)
order to control the diverging
Denoise Before Train. If the input dataset is noisy, we recommend deploying a denoising method before considering the PDE discovery task. this work, we considerod SVD and drop out components with singular values below threshold of O(10-1) to smooth out the noise. Clearly, this problem dependent. Since the noise to signal ratio is not known in advance, the commonly used denoising approaches introduce bias in the discovered PDE, regardless of deployed PDE discovery method.
29

---

Published in Transactions on Machine Learning Research (09/202))
Start
Start
Initialize parameters
Initialize parameters
0
Set
0
Find (1++1) solving system given f e
Find ((++1) by solving system PDBs given r
Find A by solving the backward adjuint eqn.
Find solving the backward adjaint
4 k
$++
Compute Gradisnt 80/ac
Update average of gradiont 8C
Update
No
Yes
(++1) T
Update parameters a
Yes
No
No
Converged?
Converged?
Yes
End
End
Figure 18: Training flowchart of the Adjoint method in finding PDEs (left) without and (right) with gradient averaging.
30

---

Published in Transactions on Machine Learning Research (09/2025)
F Impact of hyperparameters on adjoint method
In this section, we study the impact of some of the hyperparameters used in the adjoint algorithm. We repeat the PDE discovery experiment for Burgers's 's and Kuramoto Sivashinsky equation with data on a grid with N₂ = N₂ = 100, as described in 3.1.2 and 3.1.3 We check the error in the outcome coefficients and the solution of estimated forward model compared to the data.
As shown in Fig. 19 by increasing the regularization factor co, the optimization problem seems converge faster to stationary solution. In case of Burgers' equation, we considered 60 (10-1, 1000,,10=, 10-12 10-16), where for all values of 6o the exact solution is recovered. However, in the case of Kuramoto Sivashinsky with E {10-10, 100p=",,,11111 10.11110000 10010-,,- 14 10 10-16), the solution seems to be more sensitive to cp. Here, we fix the other hyperparameters Thhr = 10-16 and B = 2 x 10-3 for the Burgers's equation and 3 = 20 for Kuramoto Sivashinsky equation. We observe that high regularization factor detoriorates the accuracy, while stabilizing the regression problem.
Next, we investigate how the error changes with the thresholding tolerance where The E 110-4.10 10-8, 1100-110-11,,, 10-16). Here, fix the other hyperparameters co = 10-16 and B = 10-3 for the Burgers's equation and 3 = 20 for Kuramoto Sivashinsky equation. Although using smaller 7thrr allows faster convergence to stationary solution almost in all cases, we remind the reader that Tthr should be large enough to allow enough training of the coefficients before truncating terms. In other words, the user should avoid trivial scenarios where the initial guesses for coefficients a are zero and the thresholding is applied from the very beginning of the training.
1-11211
18 16
14
16-1
10-1
20-4
10
10
10-2
13
10-1-
50
360
320
22
75
50
Speets
Iaecta
Cooche
Opochs
(a)
(b)
(c)
(d)
1221--12
12-10-12
10-44
10-1
20
33-1
23-2
10-1
30-1
2014
900000
200
210
200
250
200
202
120
200
Specta
Iacota
S00079
Spocts
(e)
(f)
(g)
(h)
Figure 19: Impact of regularization factor 6p and thresholding tolerance Tthr on the error of adjoint method for Burgers' equation (a-d) and Kuramoto Sivashinsky equation (e-h).
Finally, we show the impact of the free parameter 3 in the learning rate on the resulting PDE discovered by the adjoint method. We compared the solution of adjoint method using 3 E (10-3,2) x 10-3, 3 x 10-3, x 10-3) for the Burgers' equation, and 3 E {2, 5, 10, 20} for the Kuramoto Sivashinsky equation. Also, we fix the other hyperparameters 7th = 6o = 10 16 As shown in Fig. 20 regardless of the value of B, adjoint method delivers the same solution. However, larger values of 3 lead to faster convergence to the solution, if the nunerical solver does not become unstable. The upper bound of 3 is limited by the stability of the guessed PDE, and can be found with try-and-error.
31

---

Published in Transactions on Machine Learning Research (09/2025)
2-26-22
5-58-02 00
6-55+00
3-20-0:
1-15
3-17
3-11
2
2-
13-1
10-4
13-1
10-1
20
00.7
20-13
10-
20-44
10-1
20
10:14
1500
(a)
(b)
(c)
(d)
Figure 20: Impact of the free parameter 3 in the learning rate the error of adjoint method for Burgers' equation (a-b) and Kuramoto Sivashinsky equation (c-d).
G Illustration of deployed notation for the considered cases.
Although the proposed method and its algorithm be and has been computed in automated fashion, here we show two detailed illustrative examples for 1-dimensional and 2-dimensional cases presented in Section 3 for the sake of better understanding the used notation and how the library of candidate terms looks like.
G.1 Heat and Burgers' Equations
As mentioned in Sections 3.1.1 and 3.1.2 for these two cases, we consider a system consisting of a single PDE, i.e. N = dim(f f = dim
where f = and = p. in a one-dimensional input space, i.e. where æ = x, and = d. In addition, we consider candidate terms consisting of derivatives with indices E (1, 2, 3) and polynomials with indices p E (1,2,3)) In other words, dmax and Pmax 3. The resulting forward model in Eq. takes the form
n = dim (x) = (d
0f C[
(AM))
Odlyy
Øt
Ord
p=1
0f
0f
Of2)
0(f*)
+ 01,1
112
01,3
0t
dr
dr
dr
8
8 f
82 (fa))
+ 02,1 0r
+ 02,
+ 0/2,3
dx2
du2
0
8 (
83 (1(
+ 03,1 0
03,2
03.3
(67)
dr
Oz
where Od,p denotes the parameter corresponding to the term with d-th derivative and p-th polynomial order. As observe, have terms with unknown coefficients (1.23)))) that aim to find using the proposed adjoint method.
The cost functional in this case is simply
C= Σ -(2)))())))) + (I,f colla j.k
(68)
Letting variational dorivatives of with respect to f to be zero, and using integration by parts, the corresponding adjoint equation can be obtained as
a
3
3
8(fP) 0d
ΣΣ
1) adP
0t
of Ord
d=1p=1
32

---

Published in Transactions on Machine Learning Research (09/2025)
0%
a1
(2
(3f22
0r
dr
dr
0º
0¹
02A
+ 02,1 Or
#222 (2
022,3 (3/2)
dx
Ox?
0³
03,1 0r
03,2 (2ʃ)
dr
a3,2
2x
(69)
with final condition {(x(0)x0+11) 2 (x(()) - f(x(0,, x(i+1))) for all j. The parameters a are then found using the gradient descent. method with update rule
ac
ad,y + ad.p Doxd.p
(70)
where 1= 3 min (A))d and ac dad.p
= (-1)² /
fo 0d Ozd
drdt + 2200dd..
(71)
This leads to the update rule for each coefficient, for example
3
0x
01, - 01,1
dxdt 2,
min
3
1² 0x
01,2 0122
min Ar)
dedl 2, dr
3
01,3 f- 01,3
dxdt 28
min Az)
dr
G.2 Reaction Diffusion System of Equations
As mentioned in Section 3.1.5 for this case, we consider system consisting of two PDEs, i.e. N dim(f) dim(p) = where f = [2] and = P2], in two-dimensional input space, i.e. n = dim(w) = dim d) = 2 where x = [2] and d = [d1, d2). In addition, we consider candidate terms with derivatives such that E Da = {[o.0) 10, [0,2]] and polynomials such that p E Dp = {[1,0] [0, 1). [1, 1), [2, 22.0).
(0,0000000000 (0, 2). [2, 1), (1, (0,3]]. In other words, dmax and Pmax = The resulting
forward model in Eq. takes the form
Liff - 8efi + Σ aidd V [fP]
(72)
d.p
where i E {1,2], fo = f4' fo and V
v ( (d2)
This is equivalent
0f:
a ta (f f1²)
Σ
+
0t
dx 0xy
[di da (pI PEIED,
0f:
+ ac.00.0).1...ff + 00.0000 + 03.10.0011.1fffff2 + ... a..0000
Øt
0h [..11.0)111.0 ox
0fz + 0..11000.00011 dr
(fiffe)
a(
+ t..(,o)
+ +... 0.. dr
dz
+
f2 8 + 0..00.2,.0011 dx² + [a. dx3
8 ( 0.,00,)),[11]
82 (
a 01.(0.2),00.3)
(73)
Oz
Oxh
where E (1,2). As we can observe, we have [Pal x P 5x 45 terms with unknown coefficients a; [ci. DD.pED for the i-th PDE, i.e. total of 90 terms for the considered system, that we aim to find using the proposed adjoint method.
33

---

Published in Transactions on Machine Learning Research (09/2025)
The cost functional in this case is simply
C=
Σ
(f) ((+11)
fi(x)(h) ((6+))))22 + M1lx, t) f)) t))dxdd colla 12.
(74)
j.
The corresponding adjoint equation is given by
04,
Σ
11
dn
PM
01
d.p
a(S9' 8 aditda d
Σ
Σ
(-1)di+ds dal 0fi
(75)
Du 0n d2
(dr |riPPaEDD
and (x() (0+1)) 2(f (x), (0+ ((( )(++1)) for all j, and where i {1,2}.
Assume, without loss of generality, that = Then, write
Qh1 (1.00.0))1..0);11 + 11.10.0).1.11] f2d1 + +a a1. 22, (20 +... + a1.00.0..3.0 (3f))
Ot
0A1 - 01,10,0,1111/1// dz 01.11,0).11.0 dz
011
dd1 a1,11.111330 (3f2)
ah1 dr
01,11,0,122000 (2f))
de
82 82 (0.01111
82 X1 +... 01.0021110 (3f)) dz3
S2 (76)
[1.[0.2).11000 Oz3
01..0.20
Oz3
dx
and (x)) (((x()))(++))) = (x ((+1) (,+(),,++)))) (*) ( for all j.k. We can follow the same procedure for = 2. The parameters ag are then found using the gradient descent method with update rule
ard,, + as.dp 7 dou,d.p
(77)
where
(-1) / ( [ccdddd + 200ar,,,p 1 min( Aa) dm- and ac
(78)
do..d.p
with Ar Ar, Ary, leading to the update rule for each coefficient, for example
04.[0.0].1.o0 + ac.]o.o).I1.o]
fid,dadt 25
min(Aa)²
9
ad, dadi - 290000,11.0000000
f1
min(Ax)
Bzi
H Setup for WSINDY
In this section, we detail the parameters used in the discovery of PDEs using WSINDy method in this manuscript following default values from Messenger and Bortz (2021) | as the reference for the reader.
Table 12: The configuration of WSINDy hyperparameters used in this work
PDE
max. derivative
max. power
(m4, m.)
(St, 8x
Burgers' Eq.
(60, 60)
(5, 5)
Heat Eq.
(60, 60)
(5, 5)
KS
(22, 23)
(5.
Random Walk
3
(60, 60)
(5.5)
34


Reviews:
[{"review_id": "Cbqup0RmN0", "reviewer": "Reviewer_dimj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is well-structured and clearly formulates its motivation to move beyond point-wise regression for PDE discovery. For readers with a background in optimal control or PDE-constrained optimization, the use of the adjoint method is an intuitive approach. However, the writing is dense and assumes significant expertise in variational calculus, which may limit its accessibility to a broader machine learning audience. While the appendices provide important derivations, the main text could benefit from a clearer bridge between the continuous mathematical theory and the practical challenges of its discretized numerical implementation, such as the stability of the chosen solvers. I did not closely check all derivations on their correctness.\n\nIn terms of method, the paper's primary contribution is the application of the established adjoint method to the inverse problem of discovering the PDE's form itself. This reframing yields an analytical gradient for the optimization, which is the source of the method's main advantages, namely that its computational cost scales with the complexity of the solver, not the number of data points, and it can naturally handle temporally sparse observations. This is in contrast to methods like PDE-FIND. However, the method's main weakness is its implementation complexity, as it requires a numerical PDE solver for both the forward and backward (adjoint) equations. Furthermore, the optimization process can be fragile, as intermediate coefficient values may define mathematically ill-posed or numerically unstable PDEs, causing the solver to fail and making the training highly sensitive to hyperparameters like the learning rate. A closer discussion of these aspects would be appreciated in the paper.\n\nThe experiments of the paper include a comprehensive suite of tests across a variety of canonical PDEs. The results suggest a high accuracy of the method and beneficial scaling on large datasets compared to PDE-FIND. The experiments also include stress tests involving noisy data, temporal sparsity, and incomplete candidate libraries. Nevertheless, the experiments are confined to synthetic data generated from known equations. Real world experimental data, which often contains complex, non-Gaussian noise, systematic biases, and uncertain boundary conditions, are not sufficiently covered. Additionally, while crucial hyperparameters are explored in the appendix, their sensitivity, particularly concerning the stability of the numerical solver, is not fully addressed in the main text, and the comparison to more recent baselines like WeakSINDy is relatively brief.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a method for identifying governing PDEs from data using an adjoint-based, PDE-constrained optimization framework. The objective is to find the coefficients within a predefined library of candidate PDE terms that minimize the L2-norm distance between the numerical solution of the parameterized PDE and the provided data. By applying variational calculus, the authors derive adjoint equations with analytical gradients of the objective function, allowing for gradient-based optimization. The reported results indicate the method can recover PDE coefficients with high accuracy and demonstrates more favorable computational scaling on large datasets compared to the baselines, though it relies on the implementation of numerical PDE solvers.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Overall the paper is already in a good form, but there are some aspects that could be improved. Firstly, the paper can be more direct about a major practical challenge, namely that the method can fail during training if it temporarily guesses an unstable equation. The authors need to explain their strategy for handling this and provide clearer, more practical guidance on how to choose the important algorithm settings and prepare noisy data. This is essential for practitioners or future research to apply the method on new data. \n\nTo further strengthen the paper, I recommend demonstrating its performance on at least one real-world experimental dataset, beyond the synthetic ones in the current paper version. It would also be beneficial to briefly explain why the method is more accurate than others in certain situations and to acknowledge the current limitation of assuming known boundary conditions. \n\nIn terms of clarifications of writing, it would be great to better explain the practical implementation of the adjoint final condition (Eq. 7). In the current paper, it is a bit unclear how this would be done in partice. Further, the precise trigger for when the algorithm starts simplifying the equation would be good to clarify.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "K0s2VAXoNg", "reviewer": "Reviewer_yG6F", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is well-written, self-contained, and mathematically rigorous. The derivation of the adjoint equations and their role in gradient-based optimization is clearly explained, and the exposition is accessible to readers outside the numerical PDE community.\n- The authors test their approach on a wide range of PDEs, including the Heat equation, Burgers’ equation, Kuramoto–Sivashinsky equation, and a 2D reaction-diffusion system, etc. This breadth of coverage enhances confidence in the method’s generality.\n- The manuscript compares the proposed method with several existing techniques, including PDE-FIND, PDE-LEARN, and WSINDy. \n- Leveraging numerical solver and adjoint methods enables exact gradient computation with respect to PDE parameters and accurate solution of the PDE.\n\nWeaknesses:\n- Although the authors test their method under additive noise, the magnitude of noise considered remains very small (typically under 0.1%). \n- All tested PDEs are relatively low-dimensional.\n- The formulation of PDE discovery as sparse regression is well established, and the adjoint method for PDE-constrained optimization is classical. The novelty of combining the two is limited—though the specific integration appears underexplored in prior work.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper cast the PDE discovery problem as a PDE-constrained optimization problem with regularization to promote sparsity in the recovered equations. The adjoint method is used to compute gradients of the loss with respect to PDE coefficients.  The resulting method allows for accurate identification of PDE structure and coefficients, and is benchmarked against several known techniques such as PDE-FIND, on many PDE problems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I have no specific revision requests at this time", "broader_impact_concerns": "I do not identify any broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "npfCpcr4oL", "reviewer": "Reviewer_8bQS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n\n* The proposed adjoint method is clearly formulated with detailed derivations.\n* The experiments are extensive, covering a lot of relevant PDE systems. However, more detailed baseline comparisons and evaluation metrics could further improve the manuscript, as I will outline below.\n\n### Weaknesses\n\nMy main concerns are about the accuracy of the experiment results, especially those related to baseline comparisons. There might be some errors in implementing the baselines, or more likely, these might be my misunderstandings due to the insufficient descriptions of experimental setups. In any case, the authors are encouraged to include more detailed descriptions and results in the experiment sections. More specifically:\n* The results for baseline comparison in Table 2 are a little surprising to me, especially the failure of WeakSINDy. Did WeakSINDy have the same hypothesis space of equations as your method? Also, could WeakSINDy do better under a different configuration from what is currently used, e.g., smaller maximum power of f and maximal order of derivatives, larger number of test functions, etc?\n* The results in Table 3 only compare against one baseline, PDE-FIND. The authors seem to justify this from the fact that PDE-FIND performs best among the baselines on the dataset of the 1D heat equation. However, since Table 3 targets a different equation, that empirical evidence may not transfer to the experiments in Table 3. Thus, I suggest the authors include other baselines as well for the experiments in Table 3. This also applies to some other tasks in the following subsections, especially: Table 9 where WeakSINDy could perform better than PDE-FIND because the observations are sparser in time and the derivative estimation by finite differences becomes inaccurate; also, the entire Section 3.3, where WeakSINDy could perform better with noisy data.\n* In the experiments for Burgers equation (Sec 3.1.2), again, the poor baseline performance is a little surprising to me. For example, the original WeakSINDy paper [1] also showed results on the 1D Burgers equation (with a different coefficient) in their Figure 4, where the true positive ratios were almost 100%. What accounts for its failure in your experiment? Is there any difference in setup?\n* Currently, each experiment on a specific dataset is run once, and the discovered equation is reported. However, the result might have some degree of randomness due to the random noise in data or the randomness in the optimization process. In many related papers, each experiment is run multiple times under different random instantiations, and additional metrics such as true positive ratio (TPR) [1], success probability, and parameter estimation error [2,3] are reported. The authors are encouraged to include these metrics to improve the statistical significance of the results when applicable.\n\nThe notations used in the paper are somewhat confusing at certain points. For example, using superscripts for powers and subscripts for partial derivatives can cause ambiguity--does $f^2_x$ mean $(f^2)_x = 2ff_x$ or $(f_x)^2$? Maybe explicitly writing out the partial derivatives (e.g. $\\partial_x f$) would be better.\n\n\n### References\n[1] Messenger, Daniel A., and David M. Bortz. \"Weak SINDy for partial differential equations.\" Journal of Computational Physics 443 (2021): 110525.\n\n[2] Kacprzyk, Krzysztof, Zhaozhi Qian, and Mihaela van der Schaar. \"D-CIPHER: discovery of closed-form partial differential equations.\" Advances in Neural Information Processing Systems 36 (2023): 27609-27644.\n\n[3] Yang, Jianke, et al. \"Discovering Symbolic Differential Equations with Symmetry Invariants.\" arXiv preprint arXiv:2505.12083 (2025).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The manuscript targets the problem of discovering governing PDEs of dynamical systems from observed solutions. It proposes to formulate the discovery problem as minimizing the difference between future timesteps in the data and predictions based on the parameterized PDE, with an additional constraint coming from the PDE form. The authors conducted extensive experiments on various dynamical system datasets, such as heat equation, Burgers equation, reaction-diffusion equation, to compare the proposed method against PDE discovery baselines such as PDE-FIND and WeakSINDy.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "* Include a derivation from eq. (3) to (4).\n* Add WeakSINDy results and details about experiment setups, as outlined in the Weakness section.\n* Add evaluation metrics, such as true positive ratio, success probability, parameter estimation error over multiple random runs, where applicable.\n* Consider changing the notations for powers and partial derivatives to avoid ambiguity.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)