D-CIPHER: Discovery of Closed-form Partial Differential Equations
Authors: Krzysztof Kacprzyk
Zhaozhi Qian
Mihaela van der Schaar
Venue: Submitted to ICLR 2023 (2023)
Keywords: differential equations
symbolic regression

Abstract: Closed-form differential equations, including partial differential equations and higher-order ordinary differential equations, are one of the most important tools used by scientists to model and better understand natural phenomena. Discovering these equations directly from data is challenging because it requires modeling relationships between various derivatives that are not observed in the data (equation-data mismatch) and it involves searching across a huge space of possible equations. Current approaches make strong assumptions about the form of the equation and thus fail to discover many well-known systems. Moreover, many of them resolve the equation-data mismatch by estimating the derivatives, which makes them inadequate for noisy and infrequently sampled systems. To this end, we propose D-CIPHER, which is robust to measurement artifacts and can uncover a new and very general class of differential equations. We further design a novel optimization procedure, CoLLie, to help D-CIPHER search through this class efficiently. Finally, we demonstrate empirically that it can discover many well-known equations that are beyond the capabilities of current methods.

Under review as a conference paper at ICLR 2023
D-CIPHER: DISCOVERY OF CLOSED-FORM PARTIAL DIFFERENTIAL EQUATIONS
Anonymous authors Paper under double-blind review
ABSTRACT
Closed-form differential equations, including partial differential equations and higher-order ordinary differential equations, one of the most important tools used by scientists to model and better understand natural phenomena. Discov- ering these equations directly from data is challenging because it requires mod- eling relationships between various derivatives that are not observed in the data (equation-dale mismatch) and it involves searching across huge space of pos- sible equations. Current approaches make strong assumptions about the form of the equation and thus fail to discover many well-known phenomena. Moreover, many of them resolve the equation-dar data mismatch by estimating the derivatives, which makes them inadequate for noisy and infrequent observations. To this end, we propose D-CIPHER, which is robust to measurement artifacts and can uncover new and very general class of differential equations. We further design a novel optimization procedure, CoLLie, to help D-CIPHER search through this class effi- ciently. Finally, we demonstrate empirically that can discover many well-known equations that are beyond the capabilities of current methods.
INTRODUCTION
Scientists have been using mathematical equations describe the world for centuries. In particular, closed-form differential equations turned out to be one of the best tools to model physical phenom- ena. differential equation describes relationship between a quantity and its derivatives (rates of change); it is called closed-form if this relationship is described by mathematical expression consisting of a finite number of variables, constants, arithmetic operations, and some well-known functions (e.g., exponent, logarithm, trigonometric functions Closed-form differential equations provide general description of reality in a concise representation that amenable to closer inspec- tion by scientists. This renders them transparent and interpretable to human experts.
Discoveries of these equations required a thorough knowledge of the theory, strong mathematical skills, substantial creativity, and good intuition. The goal of this work is to discover closed-form differential equations directly from data thus accelerating the process of scientific discovery.
Challenges in discovering differential equations from data
Partial and higher-order derivatives. Many algorithms (Brunton et al. 2016 Qian et al. 2022) only identify Ordinary Differential Equations (ODEs) which evolve only with respect one variable (usually time). In contrast, many natural phenomena are described by equations involving many variables (e.g., spatial coordinates) called Partial Differential Equations (PDEs). Many equations also involve higher-order derivatives.
Derivatives not observed. Discovering differential equations from data is challenging because the derivatives are usually not observed in the dataset (equation-daie mismatch (Qian et al. 2022 ) This makes verifying candidate equation non-trivial task. Most of the methods proposed in the literature try to resolve this issue by estimating the derivatives (Brunton et al., 2016] Rudy et al., 2017). However, estimating the derivative is difficult, especially when the data sampled infrequently or with high noise (Qian et al.) 2022 Messenger & Bortz] (021a).
Strong assumptions and constrained search space. The majority of algorithms for identifying differential equations make many assumptions about the form of the equation. In particular, they
Detailed discussion in Appendix A.2

---

Under review as a conference paper at ICLR 2023
make the evolution assumption (defined and explained later) and assume that the equation can be represented as a linear combination of prespecified functions and differential operators (Brunton ct al. 2016 Messenger & Bortz 202 However, many well-known equations, such as forced harmonic oscillator or an inhomogeneous wave equation, cannot be represented in that way.
Currently, few algorithms tackle only some of these challenges. In partieular, Weak SINDY (Mes- senger & Bortz] 2021a) is able to discover PDEs without estimating the derivative by utilizing a variational approach. However, the form of the equation is constrained be in a form amenable for a sparse regression algorithm. D-CODE (Qian et al. 2022), on the other hand, uses a varia- tional approach in conjunction with symbolic regression algorithm to discover closed-form ODEs. However, it cannot handle higher-order derivatives or multiple independent variables, so it cannot be used to discover closed-form form PDEs. The algorithms that do not require the evolution assumption appeared in (Mangan et al. 2016) and (Kaheman et al. 2020) but they require derivative estimation and only consider equations represented as linear combinations of prespecified functions.
Contributions. In this work, we develop the Discovery of Closed-form Partial and Higher-order Differential Equations a Robust Way framework (D-CIPHER) that does not estimate the deriva- tives, requires fewer assumptions, has broader search space than previous techniques, and works for both higher-order ODEs and PDEs. Our contributions are as follows:
We examine the landscape of different types of PDEs from the discovery perspective. In particular, we introduce new notions such as evolution form, evolution assumption, derivative part, and derivative-free part. We use them to describe what kinds of PDEs can be discovered with current methods and to motivate our new class of differential equations. (Section
We propose new general class of PDEs VtariationalRRuaae PDEs) that admit the variational formulation (and thus allows to circumvent the derivative estimation). We also prove a theorem that motivates a novel objective function. (Section
We use the novel objective function to develop D-CIPHER, new algorithm that searches over the Variational- Ready PDFs. (Section
We develop new optimization procedure (CoLLie) to efficiently solve a constrained least-squares problem and thus help D-CIPHER search through this space efficiently. (Section ())
PARTIAL DIFFERENTIAL EQUATIONS
In this section, we provide background information about Partial Differential Equations and intro- duce new notions necessary for the following discussion.
Notation and definitions. We denote the set 11,2....,m n as [n] and the set of non-negative integers as No. Throughout this paper we let M, N, K E be some natural numbers and let Q C Ra be an open set inside R-W. A comprehensive table with all symbols used in this work can be found in Appendix together with some definitions restated more formally.
Going beyond ODEs. The simplest differential equations are ordinary differential equations that describe quantities that evolve with respect to only one independent variable, usually time. Most methods assume that the ODE is explicit and such can be represented system of first-order ODEs:
(t = (t, unt))
(1)
where ûj represents the derivative of uj. Then the discovery problem is reduced to deciding the order of the derivative (usually first or second) and the discovery of fj.
For PDEs, it is not enough to talk about the derivative, as we can take derivatives with respect to different variables. We denote the mixed derivative as 0a where a NM is called a multi-index, and define as 8ex = 00000'> ... ...0 Each 0. = go: /0x is a a -order partial derivative with
respect to z; (the 7th independent variable | de non-trivial if [x| > 0.
We define the order of a as |ox| = EM 0g We call
A PDE of order K is any equation of the form
fax u(x), al um))) Va E Ω
Throughout this work we assume that the functions we are smooth enough for the equality of mixed partials (Spivak 2018) to hold. In that case, any mixed derivative can be uniquely specified by multi-index.
2

---

Under review as a conference paper at ICLR 2023
where @lK| u are all non-trivial mixed derivatives of all uy (j E [NT) up to the K order. We call PDE flosed-form if f is closed-form.
As PDEs might include many different combinations of derivatives, there is generally accepted counterpart of explicit ODEs in the space of PDEs.
Evolution Assumption. Although there is no generally accepted notion of an explicit PDE, define an evolution form of a PDE to be an equation of the form
(a()) = fax u(x), (lK)///(Kr)))) u(x)) Vx E E
(3)
where OlK]/a is [lk] with d omitted, ar is a known multi-index and j E [M]. Note that if M and |a]] K then Equation becomes exactly the definition of an explicit ODE.
In fact, many algorithms for PDE discovery assume particular evolution form (Messenger & Bortz, 2021a). We call an evolution assumption (EA). However, this assumption requires the knowledge of a and which might not be trivial. Usually, 8 is assumed to be the first derivative with respect time (0,) (Rudy al., 2017) but it is not the case for many well-known PDEs such the wave equation or Gauss's s law.
D-CIPHER does not need the evolution assumption. Moreover, can discover some PDEs that cannot be put into the evolution form.
Linear combinations. Current PDE discovery algorithms (Rudy et al. 2017: Messenger & Bortz, 2021a] Chen et (021)) consider PDEs that can be represented as a linear combination of functions. That means the PDE has the form
Σ Oppp(e u(x), gl u(x)) Va E Ω
(4)
where Op E IR. for porp& E [P) are the only constants that are optimized. As there are lot of expressions that cannot be put in that form, these algorithms fail to discover more complex equations. In par- ticular, for unknown 0 IR functions such sin(Or;), ex x440 cannot be learned by these algorithms.
Unlike previous methods, D-CIPHER not limited to PDEs that can be represented as a linear combination of functions. To describe the family of PDEs that D-CIPHER can discover we need to introduce the notions of derivative-bound and derivative free parts.
Derivative-bound part and derivative-fres part. Any PDE can be expressed in the form
f(x, u(x),
(l u(x)) g(æ, u(x)) 0 Væ S
(5)
where we collect all the terms with the derivatives into f(x u(x), ((x))(())) and all terms without the derivatives into g( u(x)). We call f the setivaatve-bounda part and g the derivative-fre part. We also denote them O-bound and Ò-free. The significance of the O-free part is that it can be evaluated directly given u, whereas the d-bound part requires access to the derivatives which are not observed. It is important note that for first-order ODEs, f is trivial and is just equal to úi.
The challenge of unobserved derivatives might put constraints on the d-bound part of the PDE and thus in practice, it might not be possible to search over all closed-form PDEs. However, observe that no such constraints need to be put on the Ò-free part (as it does not include the derivatives). We take full advantage of this observation and search separately for the two parts. We search over all dlosed-form functions g and for each candidate, we try to find the best counterpart f among the allowed expressions. This is very different from the previous approaches, which either do not need to find f as they work only for first-order ODEs (Qian et al., (2022) or they constrain equally both the 0-bound part and O-free part to be a linear combination of some pre-specified functions (Messenger & Bortz 2021a).
3 RELATED WORKS
Symbolic Regression. The goal of symbolic regression is to find closed-form expression that best models the given dataset both in terms of accuracy and simplicity. In contrast with the conventional regression analysis which optimizes the parameters of a pre-specified model, symbolic regression
3

---

Under review as a conference paper at ICLR 2023
Table 1: Columns correspond to challenges outlined in Section 1 and answer the following ques-
tions:
Can it discover PDEs? Does it avoid derivative estimation? Is the evolution assumption unnecessary (Equation 3 Can it discover any closed-form Ô-free part (Equation 5) References:
[1]
Brunton et al. 6) [2] (Mangan et al. 2016), [3] (Rudy et al., 2017). [4] (Long et al. 2019)
[5] Messenger & Bortz aa). [6] (Reinbold et al. 2020). [7] (Qian al. 2022)
Method
PDEs
No ô estimation
No EA
Any closed-form Ò-free part
SINDy [1]
x
x
x
X
SINDy-implicit [2]
x
X
PDE-FIND [3]
/
X
X
PDE-Net 2.0 [4]
/
X
X
WSINDy [5.6]
/
X
D-CODE [7]
X
/
X
D-CIPHER
aims at discovering both the general structure and the parameters of the model. Most of the existing work focuses on developing optimization algorithms. Genetic Programming (Koza 1992) has been widely used for that task (Schmidt & Lipson 2009). A different strategy has been employed in AI Feynman (Udrescu al. [2220] Udrescu & Tegmark 2020) that uses neural networks to reduce the search space by identifying simplifying properties like symmetry or separability. Optimization methods based pre-trained neural networks (Biggio et al. (02T), reinforcement learning (Petersen
iet al. 2021), and Meijer-G functions (Alaa & van der Schaar 2019) have also been proposed.
Data-driven discovery of closed-form differential equations. Data-driven discovery of physi- cal laws is an established area of machine learning (Bongard & Lipson 2007 Schmidt & Lipson 2009). The pioneering work in that area was SINDy (Brunton et al. 2016) that constrained the space of equations to linear combinations of functions from a predefined library and used sparse regres- sion to discover explicit ODEs. It was later extended include implicit ODEs (Mangan al., 2016 Kaheman et al.) 2020) and PDEs (Rudy et al., 2017; Schaeffer] 2017). Various other extensions were proposed by improving the derivative estimation and the training procedure (Rao et al.) 2022
Xu 2021), adding additional selection criteria (Mangan et al. (2017) and learning the library using genetic programming Maslyaev et al. 2019); Chen et al. (202I); Xu et al. (2020). A different ap- extension of (Long et al.) (018)) which uses convolutional and symbolic neural networks. is important to note that all of these methods still assume the PDE to be a linear combination as discussed in Section 2 (Equation 4 which significantly limits their search space. Some other developments are based on Gaussian processes (Raissi et al. 2018 Raissi & Karniadakis 2018) but they require the exact form of the PDE and only optimize the parameters.
proach is taken by (Long et al. 2019)
Variational approach. Recently, the variational approach has been used as a viable alternative to derivative estimation. However, they have only been used for differential equations in linear combination form (Messenger & Bortz 2021a Reinbold et al., (0200) or closed-form first-order ODEs Qian et al. 2022). Extending the variational approach to closed-form PDEs is not trivial as PDES are much more complex than ODEs and not all closed-form PDEs admit the variational formulation. In fact, the approaches that learn the library mentioned in the previous paragraph can produce exactly such terms which prohibits the use of variational formulation. To address these challenges we use the new notions defined in Section to define a new and general class of PDEs in Section 4 that admit the variational formulation.
4 VARIATIONAL-REANT
PDES
In this section, we propose a new and very general class of PDEs, the Variational- Ready PDEs (VR- PDEs), which can be characterized without referring the derivative. The VR-PDEs allow arbitrary --free part but make some minor restrictions on the Ô-bound part. These restrictions allow one to use the variational formulation of PDEs circumvent derivative estimation entirely. Despite the minor restriction, VR-PDES contain many well-known PDEs, including all linear PDEs, Maxwell's equations, and Navier-Stokes equations (additional examples provided in Appendix B.
To define the new class of PDEs, we need the following definition.
4

---

Under review as a conference paper at ICLR 2023
Definition 1 (Extended derivative and differential operator), Let No. |a| I K. be multi- index. Let : RM+N - and a RM + IR be smooth functions. An extended derivative &, denoted (ax, a, h), maps vector field R + R- to function Elu] RM + R defined as:
(() axx) [h( u(x))]
E is called closed-form if and are closed-form. We call & non-degenerate if |a] 0. Now, let (Ep)pE(e) be finite sequence of non-degenerate extended derivatives. The extended dif- ferential operator, denoted Elp is operator defined as:
((p) = Σ E (x)
Remark. Any linear operator L = EaE 0ada acting on uj is an extended differential operator. Definition 2 (Variational Ready PDE). Let Eiry be an extended differential operator, and let g : RM+N + IR be continuous function. We denote a Variational-Reut PDE (VR-PDE) by pair ()))) and define as:
(x)) g(x, uxx)) Va E
We extend the standard variational formulation of PDEs (Proposition [ in Appendix B) from linear PDEs to all VR-PDEs. The following definition is useful in further discussion.
Definition 3. Consider field u Q + RN,, and an extended derivative & (a,, a, h). Let
6: . + be a testing function (C function with compact support). We define the functional
F(E, 6) - h(x, ((z))(())((1)) o (a(2))(())d (((
We can use this functional to formulate variational characterization of VR-PDEs.
Theorem 1. . RN, where uj E CKK, is solution to VR-PDE in Equation if and only if
P
F(Ep, u. Φ) - -f [9(0, [loto u(=)))(x)) dx =
Σ
(9)
p=1
for all testing functions 0:02 RR
Proof. Appendix
This theorem motivates the variational loss function as we expect the left-hand side of Equation to be closer to 0 the closer the canditate PDE to the true one. To calculate how well a set vector matches VR-PDE (E11).9) propose the following loss function.
fields D (u( Id=1
D s P (()
2
C (E)p))9 = ΣΣ (l
(.)
g(2, (d) (2))
(10)
d=1s=
A
where {0.)s- is a set of predefined testing functions.
This novel loss function makes it possible evaluate to what extent any VR-PDE matches the observed data. This loss can be used as an optimization objective in any algorithm that searches over some subspace closed-form VR-PDEs. We propose D-CIPHER in Section 5 as an example of such an algorithm.
D-CIPHER
In this section, formulate the problem of PDE discovery and then introduce novel algorithm (D-CIPHER) to solve it.
We say u: R R is in CK 8 exists and is continuous for all |a]| S K.
5

---

Under review as a conference paper at ICLR 2023
Problem formulation We given a dataset of observed fields D -(b(d) with a finite sam- pling grid 9 C S2. Each w(()) (a:) is noisy measurement, i.e., v(d) g - IR N is defined
(a) (d) (d) (x) "j (xx) + + (a) Vo GVj E [N]
(1)
(11)
(d) Ω + R is C where ( (az) is realization of a zero-mean random variable (noise), each 3 function, and every true field u(d) is governed by the same closed-form PDE f. The task is infer the closed-form PDE. f. from the dataset D {(r(9) | and the sampling grid G. We assume that I is inside the class of closed-form VR-PDEs (Section 4 and its 0-bound part is inside a subspace of extended differential operators spanned by user-specified dictionary (see Step below).
We propose an algorithm that consists of three steps. In the first step, we define the subspace of closed-form VR-PDEs we want search over to reflect our knowledge of the problem. In the second step, we reconstruct the fields from noisy measurements. In the last step, we solve an optimization problem using modified symbolic regression algorithm. For more details, check Appendix C
Choose the form and incorporate prior knowledge. A human expert should encode their prior knowledge of the problem into dictionary of non-degenerate extended derivatives Q [E)] We use this dictionary to search over a finite mime--iimeniionl subspace of closed-form operators spanned by this set. In other words, we assume that the VR-PDE is of the form:
Step 1:
spp(((u)(z) g(x u(z)) Va Ω
(12)
where B E IR.", is any dlosed-form functiin of M + N variables, and E (app ap. hp).
For instance, dictionary might include only the partial derivatives up a certain order. For 1+1 second-order equation that means Q = {8,, Ôz 8zz, 00,0% j.). That is already enough to discover heat and wave equations with any closed-form source. If, for instance, the user suspects the presence of the advection term udy (as in the Burgers' equation), the term 0g (u²) can be included in the library.
It's important to note that we do not assume any particular form of g apart from being closed-form.
Step 2: Estimate the fields. As the dataset D consists of noisy and infrequently sampled fields, we first need to estimate the "true" fields (d) from b(d), Any choice of reconstruction algorithm can be used and the user should choose it according the problem setting and their domain knowledge.
Step 3:
Optimize. We minimize the loss function in Equation 10 for the estimated fields (ā( ) among all PDEs of the form in Equation 12 We solve the following optimization problem:
D
S
min
min |||||1
Σ pml
F(B a .b)))
g( (
(13)
As we want to discover both g and B we cannot use the standard penalties on B such as the \||,8||2 or 111 as the loss would be minimized by O and B O. Therefore we put the constraint ||@||1 = 1. We choose the LI norm to encourage sparsity in the coordinates of the vector 3.
The inner minimization in Equation
13 can be rewritten as constrained least-squares problem.
min ||
Σ
(3 (3.z(d.)) w(d.))
(14)
(d.s)
where E, (Crp, ap hp) and (dd,)) E IR,,, wddde) E R are defined as
2
(d,s)
hp (x, (()) (z)) (a)))))))))
(a
1
(15)
(d.s) / g(1r (((0 (())))))))) (z))
We show the full derivation in Appendix C 2 (d.s) can be precomputed at the beginning of the algo- rithm without estimating the derivatives of the reconstructed fields. They can be easily calculated the derivatives of the testing functions ob and the derivatives of a, can be analytically computed.

---

Under review as a conference paper at ICLR 2023
As the optimization problem in Equation 14 has to be solved many times for different closed- expressions g. it poses some unique challenges. As standard approaches are not sufficiently fast, we design new heuristic algorithm to solve this problem. We describe it in the next section.
COLLIE
The problem in Equation 14 from the previous section can be formulated follows. Given matrix A E Im*x" and vector b Rm find a vector & E II" that minimizes | such that ||=||| 1.
The task is challenging as the unit I.1 sphere is convex. A method that guarantees an optimal
solution is based on an observation that the (n 1)-dimensional L.1 sphere consists of 2 (n
simplices (which are convex). Minimizing A bl12 on simplex a quadratic program (Boyd
ct al. 2004) with many available solvers (Andersen et al. 2013 Stellato el al. 2020 ApS 2019).
However, that means that the computation time scales exponentially with the number of dimensions.
This is prohibitively long for the inner optimization of our algorithm. Therefore we design a heuristic algorithm CoLLie (Constrained LI Norm Least Squares) that finds an approximate solution but
significantly faster (Figure I. We provide a detailed description of CoLLie in Appendix D
(B)
(C)
100
Avg. relative error
0.009
Relative error
10"
10-3
10-1
CoLLie CVXOPT
10
10-19.
10-3
-15
5
0.0
Figure 1: We compare CoLLie with an algorithm that uses CVXOPT (Andersen et al. (013)) to solve each of the convex subproblems. We report the relative error between the loss obtained by CoI.Lie and the minimum loss achieved by CVXOPT. Panels B and C show the averages and the distributions of relative errors. The average relative error below 0.005 and the bulk of the distribution below 10 At the same time CoLLie is orders of magnitude faster (Panel A).
7 EXPERIMENTS
We perform series of synthetic experiments to show how well D-CIPHER is able to discover some well-known differential equation (Table 2). First, we demonstrate that D-CIPHER performs better than current methods when discovering PDEs in linear combination form (Section 7.1). Then demonstrate it can discover PDEs with a closed-form Ò-free part that cannot be expressed as a linear combination and thus are beyond the capabilities of current methods (Section 7.2) We contrast D- CIPHER with its ablated version where the derivatives are estimated and the standard MSE loss is used instead of the variational loss (details in Appendix E. For additional information about the experiments (e.g., implementation details, data generation, experimental settings) see Appendix
Table 2: Equations used in the experiments. "LC" column specifies if the equation can be repre- sented linear combination (Equation "VR" column specifies if the PDE Variational- Readd
Name
Equation
LC
VR
Homogencous heat equation
0,u 0 0.uu 0
1
/
Burger's 's equation
0,u + uº,u - -002u 10Pu=0 =0
/
/
Kuramoto- Kuramoto-Sivashinsky equation
0,u + Pu + -.uu -uduu=
/
Forced and damped harmonic oscillator
0f 20 0z0,u Ozu Og sin(O,t)
X
SLM model (Appendix F.I
d,u + O.
-2
X
Inhomogeneous heat equation
0 0 OL = Ost
X
Inhomogeneous wave equation
8L 01 02 u Oze' sin(Agt)
X
Evaluation metrics. To establish how well a discovered PDE matches the ground truth, we evaluate its Ò-free and 0-bound parts separately. For the Ò-free part, we assign binary variable indicating
All experiment code will he published upon acceptance.

---

Under review as a conference paper at ICLR 2023
whether the correct functional form of the equation recovered (please check Appendix E.8 for details). For the d-bound part, we measure the RMSE between the found coefficients of 3 and the target ones. We report the averages and standard deviations for both parts. We call the averages respectively Success Probability and Average RMSE..
Implementation. We use B-Splines (De Boor 1978) the testing functions and we estimate the fields in Step 2 of D-CIPHER with a Gaussian Process (Williams & Rasmussen (006). The outer
optimization in Step is performed using modified genetic programming algorithm (Koza 1992) and the inner optimization by CoLLie (Section 6)
7.1 DISCOVERING LINEAR COMBINATIONS: COMPARISON WITH OTHER METHODS
We compare D-CIPHER against two variants of PDE-FIND Rudy al. 2017 and WSINDy Rein- bold et al. 2020) with optimization performed by Stepwise Sparse Regression [Boninsegna et al, 2018) Forward Regression Orthogonal Least-Squares (Billings) 2013 We note that D-CIPHER is specifically designed to discover PDEs that are beyond the capabilities current methods, i.e., where the derivative-free part can be any closed-form expression. Current methods usually tested on equations where the derivative-free part is trivial (identically equal to 0). Even though these al- gorithms are specialized discover these simpler kinds of equations, D-CIPHER performs better than (or equally well as) PDE-FIND and WSINDy, regardless of the optimization algorithm, when tested on Burgers' equation the homogencous heat equation, and Kuramoto- Siviittoo-Svssskkk equation (Figure 2). This demonstrates gain from both the variational loss and the new optimization routine.
Homogeneeus heat equation
Burgers' equation
Surameto-Sivashine equstion
0.5
0.25
BMSB
0.08
0.1
WSiNDy (SIRIF)
Rvg RMS
0.4
Avg. RMSE
0.2
PDE AND (SI8FF)
D-CIPHER
0.3
0.15
PDE-FIND (F)
0.04
0.2
0.1
WSNDY (SW)(F)
WSINDY (SI
0.02
POE AND ISIRIF)
0.1
POE- TOE-RINDIO (S14F)
0.05
WSINDY
10
10
10
0.0
10-3
10
10
0.0
10-
10-2
10
Noise ratio
Noise ratio (oa)
Noise ratio lox)
Figure 2: Simulation results for the Burgers' equation, homogeneous heat equation, and Ku-
Sivass-Svv equation. We report the Average RMSE of the Ô-bound part of the equation. Note that some of the benchmarks overlap
7.2 DISCOVERING EQUATIONS BEYOND CURRENT METHODS
Forced and damped harmonic oscillator. As the oscillator is described by a second-order ODE, cannot be discovered by D-CODE (Qian al. 2022). D-CIPHER discovers the correct functional form of the 0-free part and achieves a low RMSE for the coefficients of 3 in most of the experimental settings. The performance is higher than comparable to the ablated version of D-CIPHER, thus
demonstrating gain from using the variational approach. We present the results in Figure 3
Forced and damped harmonic oscillator
Sucoess Probabiity
0.8
Succass Probabiity
Prebabiity
0.8
0.6-
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.0
1030
10
10a
0.3
15
D-CIMER
Noiseratio lou)
Sampling interval 1841
Number of samples (0)
0.08
0.08
0.00
Rvg, ReASE
0.06
0.08
Avq.
Avg.
0.02
0.04
0.0
30-4
IdP
0.0
Noise rabi
Sampling Interval 1ds1
Number of samples
Figure& 3: Success probability of discovering the correct O-free part the equation and the average RMSE between the recovered O-bound part and the target one across different experimental settings. We compare D-CIPHER against its ablated version (Abl. D-CIPHER).
8

---

Under review as a conference paper at ICLR 2023
Inhomogeneous heat equation. D-CIPHER is able to discover the correct equation in settings with very high noise. It performs better than the ablated version, thus showing the importance the variational objective. The result are presented in Table
Table We report the success probability of discovering the Ò-free part and the Average RMSE of the 0-bound part for the inhomogeneous heat equation. Standard deviations shown in brackets.
Success probability
Average RMSE
Method
0 0.05
0.1
0.2
0.05
0.1
0.2
D-CIPHER
0.64 (.07)
0.42 (.07)
0.12 (.05)
0.15 (.009)
0.21 (.007)
0.24 (,005)
Ablated D-CIPHER
0,46 (.07)
0.20 (.06)
0.04 (.03)
0.18 (.009)
0.24 (.008)
0.27 (.007)
Inhomogeneous wave equation. This equation does not have the standard evolution form, as does not involve the o term. Thus, even without the source term, most of the current methods cannot be applied directly to discover this equation. In Figure 4 we show the absolute difference between the true field and the fields computed from the sources discovered by D-CIPHER and its ablated version across different measurement settings. D-CIPHER finds the correct functional form with coefficients far from the ground truth. The ablated version fails to discover the correct functional form and the found Ò-free part does not reproduce the correct behavior of the equation.
10
10
2.02 elsin(2.98d
2.11 sint2.977
2.17 x e sints 94t)
Error
max
D-CIPHER
2.52 12 sinte")
2.99 x !sinte") -9
2.16 log xej]
Abl. D-CIPHER
Figure 4: We solve the inhomogeneous wave equation for the Ò-free parts found by the D-CIPHER and its ablated version (Abl. D-CIPHER). We show the absolute difference between the computed fields and the true field generated by Ò-free part et sin(3t).
8 DISCUSSION
Applications. As D-CIPHER can potentially discover any closed-form Ò-free part, especially useful when this part of the PDE captures essential component of the phenomenon. We demon- strate it by finding the heat and vibration sources as well as the driving force of an oscillator. Beyond the spatio-temporal physical equations, D-CIPHER might prove useful in discovering population models structured by age, size, and spatial position (Webb) 1985 2008), age-dependent epidemio- logical models Hoppensteadd (1974). and predator-prey models with age-structure (Promrak al.) 2017). All these equations are VR- PD-P where the Ò-free parts are crucial elements of the equations signifying the rates of mortality, infection, recovery, or growth.
Limitations and open challenges. D-CIPHER may fail in some scenarios, either due to challenging experimental settings a challenging underlying PDE. Challenging experimental settings might include unobserved variables, high measurement noise, infrequent sampling, and inadequate domain (e.g., small time horizon). Challenging PDE forms might include a PDE outside the VR-PDE class or O-free part with a complex expression that is difficult to find. We note that we address some of these challenges by utilizing a variational approach, defining VR-PDEs to be a very general class of equations, and designing CoLLie enabling thorough search across closed- form expressions.

---

Under review as a conference paper at ICLR 2023
Ethics Statement. We want to emphasize that D D-CIPHER designed to facilitate the process of scientific discovery by extracting closed-form. PDEs from data. is not intended to or capable of replacing human experts in the modeling process. No human-derived data was used.
Reproducibility Statement. The assumptions of Theorem T are discussed Section and the proof is presented in Appendix [B.2] The details of D-CIPHER and CoLLie discussed in Section 5 and as well as in Appendix and D The implementation details, data generation procedures, hyperparameters, and experimental settings are described in Appendix for D-CIPHER and in Appendix D for CoLLie. All experiment code will be published upon acceptance.
REFERENCES
Ahmed M Alaa and Mihaela van der Schaar. Demystifying black-box models with symbolic meta- models. Advances in Neural Information Processing Systems, 32:1 92:11304-11314 2019.
MS Andersen, Dahl, and L Vandenberghe. CVXОРT: A Python package for convex optimization, Version 1.1. 6, 2013, 2013.
MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual. Version 9.0., 2019. URL attp://docs mosek. com/ 9, wtt.//ooooo/ooee.......oooeeo html
L. Biggio, T. Bendinelli", A. Neitz, A. Lucchi, and G. Parascandolo. Neural Symbolic Regression that Scales. In 38th International Conference Machine Learning, July 2021.
Stephen A. Billings. Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Peatio-Temporal Domains. John Wiley & Sons, September 2013. ISBN 978-1-119-94359-4
J. Bongard and H. Lipson. Automated reverse engineering of nonlinear dynamical systems. Pro- ceedings of the National Academy of Sciences, 104(24):9943 0e009999999999,,,,0 June 2007. ISSN 0027-8424, 1091-6490. doi: 10.1073/pnas. 0609476144044
Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi. Sparse learning of stochastic dynamical equations. The Journal of Chemical Physics, 148(24):241723, June 2018. ISSN 0021-9606,
1089-7690. doi: 10.1063/1.500 85088
Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge uni- versity press, 2004.
Richard Brent. Algorithms for minimization without derivatives. Courier Corporation, 2013.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the Na- tional Academy of Sciences, 1 113(15):3932-390 -937,, April 2016. ISSN 0027-8424, 1091-6490. doi:
10.10733pnas.1551 151.///4...117381115.. 1517384113.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Yuntian Chen, Yingtao Luo, Qiang Liu, Hao Xu, and Dongxiao Zhang. Any equation is forest: Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). arXiv:2106.1 1120006.199227 [cs.NE], 2021.
D Crighton. Model Equations of Nonlinear Acoustics. Annual Review of Fluid Mechanics, 11(1):
11-33, 1979, doi: 10. 16/////rrev.f.... 11..111.0000 11.010199.000300
Carl De Boor. A practical guide to splines, volume 27. veringer-verlag New York, 1978.
Bradley Efron, Trevor Hastie, lain Johnstone, and Robert Tibshirani. Least Angle Regression. The Annals of Statistics, 32(2):93, 2004.
Lev Elsgolc. Calculus of variations. Courier Corporation, 2012.
G. (Friedrich Gerard) Friedlander. Introduction to the theory of distributions / F.G. Friedlander. Cambridge University Press, Cambridge, 1982. ISBN 0-521-24300-9.
10

---

Under review as a conference paper at ICLR 2023
Frank Hoppensteadt. An Age Dependent Epidemic Model. Journal of the Franklin Institute, 297 (5):325-333, May 1974. ISSN 00160032. doi: 10. 001101690166000322499990344444
James M. Hyman and Basil Nicolaenko. The Kuramoto-Sivashing equation: A bridge between PDE'S and dynamical systems. Phyvica D: Nonlinear Phenomena, 18(1):113-126, January 1986.
ISSN 0167-2789. doi: 10.1016/0167-278998686899
Kadierdan Kaheman, J. Nathan Kutz, and Steven L. Brunton. SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Soci- ety A: Mathematical, Physical and Engineering Sciences, 476(2242):20200779, October 2020. URL htt.appes/rrpaaoooeeyyypoiiiiigo
doi: 10. 10.10988rspa.2020. wtw.:p..............22 0279.
10 10.1099//1 rspa. 2020. 0279 Publisher: Royal Society.
Doyo Kereyu and Genanew Gofe. Convergence Rates of Finite Difference Schemes for the Diffusion Equation with Neumann Boundary Conditions. American Journal of Computational and Applied Mathematics, pp. 11, 2016.
John R. Koza. Genetic programming: on the programming computers by means of natural selec- tion. Complex adaptive systems. MIT Press, Cambridge, Mass, 1992. ISBN 978-0-262-1111 117-6-
Yoshiki Kuramoto. Instability and Turbulence Wavefronts in Reaction-Diffusion Systems. Progress of Theoretical Physics, (866):885-1999) June 1980. ISSN 0033-068X. doi: 10.1 10.1143/ PTP.63.1885.
Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai. When gaussian process meets big data: A review of scalable GPs. IEEE Transactions on Neural Networks and Learning Systems, (1 4405 4423, 2020. doi: 10.1 10.1109/TNNISS. wwt............. 299199997999
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-Net: Learning PDEs from Data. In Proceedings of the 35th International Conference on Machine Learning, pp. 3208-3216. PMLR, July 2018.
ISSN: 2640-3498.
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, 108925, December 2019. ISSN 0021-9991. doi: (0006//..pp2219.10889998
Lynn H. Loomis. Advanced calculus Lynn H. Loomis and Shlomo Sternberg. Addison-Wer Wesley series in mathematics. Addison- Addison-Wesley Pub. Co., Reading, Mass. London, 1968. ISBN 0-201- 04305-X. Publication Title: Advanced calculus.
N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor. Model selection for dynamical systems via sparse regression and information eriteria. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2204):20170009, August 2017. ISSN 1364-5021, 1471- 2946. doi:: 10. 0t998//spa.201..0000
Niall M. Mangan, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Inferring Biological Networks by Sparse Identification of Nonlinear Dynamics. IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2(1):52-63, June 2016. ISSN 2332-7804. doi: 10. 1109/TMBMC.201 2.663326.. Conference Name: IEEE Transactions on Molecular, Biological and Multi-Scale Communications.
Mikhail Maslyaev, Alexander Hvatov, and Anna Kalyuzhnaya. Data-Driven Partial Derivative Equa- tions Discovery with Evolutionary Approach. In João M. F. Rodrigues, Pedro J. S. Cardoso, Jânio Monteiro, Roberto Lam, Valeria V. Krzhizhanovskaya, Michael H. Lees, Jack J. Dongarra, and Peter M.A. Sloot (eds.), Computational Science ICCS 2019, Lecture Notes in Computer Sci- ence, pp. 635-641, Cham, 2019. Springer International Publishing. ISBN 978-3-030-22750-0. doi: 11110009778-3000-222555-
Daniel A. Messenger and David M. Bortz. Weak SINDy for partial differential equations. Journal of Computational Physics, 143:110525, October 2021a. ISSN 00219991. doi: 10....6/f..ep.222.. 110525.
11

---

Under review as a conference paper at ICLR 2023
Daniel A. Messenger and David M. Bortz. Weak SINDy: Galerkin- -asee lection. Multiscale Modeling & Simulation, 19(3):1474-144 1497, January 2021b. ISSN 1540-3459, 1540-3467. doi: 10.1137720M134316
Model Se-
Aaron Meurer, Christopher Smith, Mateusz Paprocki, Ondiej Certík, Sergey B Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, and others. SymPy: sym- bolic computing in Python. Peerl Computer Science, 3:e103, 2017. Publisher: Peer] Inc.
A. R. Mitchell and D. F. Griffiths. The finite difference method in partial differential equations. Wiley, Chichester [Eng.]: New York, 1980. ISBN 978-0-471-27641-
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-leam: Machine Learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne Kim. Deep Symbolic Regression: Recovering Mathematical Expressions From Data via Risk-seeking Policy Gradients. ICLR 2021, 2021.
Promrak, G. C. Wake, and C. Rattanakul. Predator-pre prey Model with Age Structure. The ANZIAM Journal, 59(2):555-166, October 2017. ISSN 1446-1811, 1446-8735. doi: 10.1017/ 0146188117000366 Publisher: Cambridge University Press.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-CODE: Discovering Closed- form ODEs from Observed Trajectories. The Tenth International Conference on Learning Repre- sentations, 2022.
Maziar Raissi and George Em Karniadakis, Hidden physics models: Machine learning of nonlin- partial differential equations. Journal of Computational Physics, 357:125-141, March 2018. 8. ISSN 0021-9991. doi: 10.10161).ee2221111..33
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Numerical Gaussian Processes for Time-Dependent and Nonlinear Partial Differential Equations. SIAM Journal on Scientific Computing, 40(1):A172-A198, January 2018. ISSN 1064-8275, 1095-7197. 17M1
doi: 10.1137/
Chengping Rao, Pu Ren, Yang Liu, and Hao Sun. Discovering Nonlinear PDEs from Scarce Data with Physics-Encoded Learning. ICLR 2022, 2022.
Patrick A. K. Reinbold, Daniel R. Gurevich, and Roman O. Grigoriev. Using noisy incomplete data to discover models of spatiotemporal dynamics. Physical Review E, 101(1):010033, January 2020. ISSN 2470-0045, 2470-0053.
doi: 10. 10.1100/yyyse 101.010203..
Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4):e1602)) April 2017. ISSN 2375-2548. doi: 10. 11111//ccady.16660 1602614.
Hayden Schaeffer. Learning partial differential equations via data discovery and sparse optimiza- tion. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473 (2197):201 60446, January 2017. ISSN 1364-5021, 1471-2946.
doi: 10. 10. (0.1098/rspa.201 2.06444
Michael Schmidt and Hod Lipson. Distilling Free-Form Natural Laws from Experimental Data. Science, 324(5923):81-85, April 2009. ISSN 0036-8075, 1095-9203. 1165893.
doi: 10. 20.1126/cceence.
Michael Spivak. Calculus On Manifolds: Modern Approach To Classical Theorems Of Advanced Calculus. CRC Press, 2018, ISBN 978-0-429-50190-6 077-:-019999097----------9999999 978-00-49-970044550 OCLC: 1029237047.
B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd. OSQP: operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4):637-672 2020. doi: 10. 10.10077/12532-0200000
12

---

Under review as a conference paper at ICLR 2023
Trevor Stephens. gpleam: Genetic programmingg in python, with sciki- inspired and compati- ble upi, 2022. URL httpt::/pppeaartt//// readt http:///ww.....tot//ooo/eeaaaaaa stable /indllee//edxx index. html
Robert Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical
Society. Series (Methodological), 58(1):267-288, 1996. ISSN 0035-9246. Publisher: [Royal Statistical Society, Wiley].
Silviu-Marian Udrescu and Max Tegmark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16)eaay2631, April 2020. ISSN 2375-2548. doi: 10.1
sciadv.aay2631.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. 34th Conference Neural Information Processing Systems (NeurIPS 2020), 2020.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy. David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, Jarrod Millman, Nikolay Mayorov, Andrew R. Nel- son, Eric Jones, Robert Kern, Eric Larson, C J Carey, Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mul- bregt. and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272. 2020. doi: 10.1038/s4 1592-119--68888888888
Sachin Wani and Sarita Thakar. Crank-Nicolson Type Method for Burgers Equation. International Journal of Applied Physics and Mathematics, 3:324-328, January 2013. doi: 10.7763/IAAPM
2013.V3.230.
G. F. Webb. Population Models Structured by Age, Size, and Spatial Position. In M. Morel, F. Takens, B. Teissier, Pierre Magal, and Shigui Ruan (eds.), Structured Population Models in Bi- ology and Epidemiology, volume 1936, pp. 1-49. Springer Berlin Heidelberg, Berlin, Heidelberg. Series
2008. ISBN 978- 978-3-5400-8277-88 978-3-540-78273-55 doi: 10.
Title: Lecture Notes in Mathematics.
Glenn Webb. Theory of Nonlinear Age-Dependent Population Dynamics. CRC Press, January 1985. ISBN 978-0-8247-7290-1.
Christopher K Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Hao Xu. DL-PDE: Deep-Learning Based Data-Driven Discovery of Partial Differential Equations from Discrete and Noisy Data. Communications in Computational Physics, 29(3):698-728, June 2021. ISSN 1815-2406, 1991-7120. doi: (0.4208/cicp.0A-2220-00
Hao Xu, Haibin Chang, and Dongxiao Zhang. DIGGA-PDE: Discovery of PDEs with incomplete candidate library via combination of deep learning and genetic algorithm. Journal of Computa- tional Physics, 318:109584 October 2020. ISSN 0021 1-991.
doi: j0.1016/j..cp.20200 109584.
13

---

Under review as a conference paper at ICLR 2023
TABLE OF SUPPLEMENTARY MATERIALS
1. Appendix
notation and definitions
2. Appendix
variational formulation for linear PDEs and the proof of Theorem
||
3. Appendix details of the D-CIPHER framework, including pseudocode
4. Appendix D details of the CoLLie algorithm
5. Appendix details of experiments and the implementation
6. Appendix additional experiments and discussion
A NOTATION AND DEFINITIONS
A.1 NOTATION
Table 4: Symbols used in this work
Symbol
Meaning
[A]
set of numbers (1....,
N
set of natural numbers, i.c. {1, 2, .....
No
a set of non-neee non-negative integers, ic.. 1, ..... the dimension of the domain of vector field the dimension of the codomain of a vector field denotes the smoothness of functions the maximum order of derivatives the size of the dataset of observed fields
M
N
K
D
S
the number of testing functions
Ω
open set in IR the derivative of u at multi-index, an element of N&
i(t)
a
la
the order a. ax = Σ M a;
00:
th -odder partial derivative with respect the 7th variable 8 022 ... 0
g
ck
a set of functions with continuous partial derivatives d for all aa| S K an extended derivative, Definition |
Elp
an extended differential operator, Definition 1
F
the functional used in the variational loss, Dehnition
(()
the variational loss, Equation 10
g
sampling grid, Definition 10
true field, Definition an observed field, Definition 10
v
D
dataset of observed trajectories the noise
Q
dictionary of non-degenerate extended derivatives vector describing the d-bound part of the VR-PDE
3
OR
noise ratio
At
sampling interval
A.2 DEFINITIONS
In this section, we collect the definitions of some of the important terms used in the paper for easy reference.
Definition (Closed-form expressions and functions). A closed-form expression is mathematical expression that consists of finite number of variables, constants, arithmetic operations, and certain well-known functions (c.g., logarithm, trigonometric functions). A function f called flosed-form if it can be represented by closed-form expression. F.g.. far y) x² log(y) + sin(3z).
14

---

Under review as a conference paper at ICLR 2023
Remark. In practice, do not want to consider any finite expression. Any symbolic regression algorithm penalizes expressions that are too long putting a soft constraint on the number of elements used. That is why deep neural networks are not considered closed-form even if they satisfy the conditions in Definition
Definition (Mulli-index).. An n-dimensional muulti-index a is an n-tuple
ax = (0r1, 022,- 02,-..,,o. On)
where Vi E [n]] ag No. Thus No. We define the order of a |a| > az.
Definition 6. For any n-dimensional multi-index a we define a mixed derivative
8 8 ".." 0
where 8h = 8a: /0x, a -order partial derivative with respect to z; (the 7th independent variable). We call da non-trivial if [x] > 0. We denote the list of all non-trivial partial derivatives of u up order K 0.k)u.
Definition (Closed-form Partial Differential Equation). Let f be closed- real smooth func- tion. We say that vector field u Q > TRN is governed by a W -order closed-form PDE described by if
fa, u(x), (((x),kkk)))) alk u ()) Va E E
(16)
where OlK) u are all non-trivial mixed derivatives of all Uy (j [NT up the T order. Definition 8 (Testing function), Support of a function Φ - R is defined as
supp Φ = {ar Q: ((x) 7 0)
where B is the topological closure of B in Ω2
Φ is called a testing function if it is C function with compact support.
Definition 9 (True Field). We define a true field on Q vector valued function u R.V where each Uy Ω > IR is a CK function.
Definition (Observed field and sampling grid). We define sampling grid G to be finite subset
of Ω. Let Ω -> R.W be true field ft An observed field sampled from grid is function v g + R of the form
Vy(x) = (x( + (z( Va E G GVj E [N]
where (x) corresponds to noise, a realisation of a zero-mean random variable.
Definition 11 (L1 sphere). Let N. We define n-dimensional L1 sphere be a subset of 1=+1 defined as:
(ar 1R*+1 ||ع||| = 1} c m=*1
(17)
Definition 12 (Standard simplex). Let N. We define standard n-simplex to be subset of IR defined as:
a+1
fx E Ra+1
I = 1 2 20EE [n 1}} C Ra+1
(18)
B VARIATIONAL-READY PDES
B.1 VARIATIONAL. FORMULATION OF PDES
In this section, we provide the standard variational formulation of PDEs for linear PDEs (Fredlaan der 1982).
Definition 13 (Linear differential operator). Let A be finite set of iultiiindiies.. A linear differ- ential operator L is defined
L =
Σ 0a aEA
where ax CK is non-zero sufficiently smooth function of dependent variables. If max |a n then we call L an th order linear differential operator. If all a are constants we say that L has constant coefficients.
15

---

Under review as a conference paper at ICLR 2023
The adjoint of L, denoted L1, is linear differential operator defined as
L ul(x)= = ' (aa))
(19)
aEA
Proposition 1 (Variational Formulation of PDES for linear PDES), Let K E N. Consider a scalar field u Ω I. such that u E CK₁, -oddr linear differential operator and continuous function S2 + R. Let S2 + R. be a testing function. Then satisfies a linear PDE
L|u(x)] - g(x) -V E
(20)
and only if
1 (u(x)[(((a) {()))(x)) dx =
(21)
for all testing functions Φ 92 R.
Note that the integrals are always well-defines as & has compact support.
B.2 THEOREM
Before we prove the Theorem [ we need the following lemma, which is a particular formulation of the Fundamental lemma of calculus of variations (Elsgole 2012). We also need generalized version of the divergence theorem (Loomis) 1968).
Lemma (Fundamental lemma of calculus of variations). Let K N, S2 be an open set in R.M. and Ω + I. be continuous function. Then u. is equal to 0 on the whole Q if and only if Sa (())(())da for all CK functions
- with compact support.
Proof. If u is identically 0 on Ω then all the integrals trivially equal to We now prove the converse.
Let assume for contradiction that there exists point xo Ω such that u((zt)) * O. Without loss of generality we assume u (2)) = 6 > 0. As u continuous there exists an open ball around Do of radius 8, denoted 1S, {x Ω || æoll2 d}. such that Væ DS u(x) e/2 >
Now let Φ be C function that is positive on B and elsewhere. Such function can always be created by appropriately shifting and scaling (xx) -11((111==11) [||=||<11] Its support is closed ball B6 (æ Ω || æol12 8) which is compact. Then
((())(()))d
u(=))(())d
(22)
n
Bé
as both u(x) and (xx) are positive on Bá een support such that:
Thus we found a continuous function & with compact
L u(z)((x)dd +
(23)
Therefore u is identically 0 on Ω.
To make this section self-contained, we provide the statement of the generalized divergence theorem (Loomis 1968).
Theorem 2 (Divergence theorem). Let Ω be open in M and let f,g be continuous = U N2 and continuously differentiable on Q. Then
[ 0 (( =- f(x)) ff v)f(z)g(z))d
(24)
where D is normal unit vector to the boundary OΩ.
In 1-dimensional setting, the statement of the theorem reduces to the integration by parts.
We prove Theorem
16

---

Under review as a conference paper at ICLR 2023
Proof. Let denote Ep (ap, ap. hp). Then the PDE in Equation 8 can be written
P
Σ
(p())) (p(a))00x1ppp1
(2 u(x)) g( u(x)) = 0 Var E
(25)
The LHS is continuous as all dp and hp are smooth, is continuous, u E C and |app S K Vp [P]. Thus we can use the fundamental lemma of calculus of variations say that the Equation 25 is true and only if
Σ (x)x)) u(x))] g(x, u(x))) Φ(x)dx
(26)
for all testing functions p. We transform the LHS of Equation 26 using linearity to:
ap xp(x)00*1hpp (ap u(x)))o(x) / g(x, u(z)))(())))
Σ
(27)
Let us now focus on
/ (p(2)00 (2p u()))((()))
(28)
and let us denote ap (apl..... (pp)....,ppp))) Then 0ap 8pm S Shg and the expression can be written as
/ 8fmh ....... .r.... [ap(a, u(x)) [ap ()))))))
(29)
Let us denote the support of @ as B. As equal to zero outside of its support, we can write the expression as
0p ...
...0%p 0.p [,p( (z))))a((z))(()))
(30)
Without loss of generality, let assume that Opl 0. By the divergence theorem, this can be rewritten as
/ 0p--- .. (, ((
(31)
because the integral over the boundary is equal to
110 0..0% (2, ul u( =
(32)
as
& has compact support (and thus vanishes the boundary). We can perform this operation dp1 times shift the whole derivative 0.n to the second part of the equation and obtain
(-1)0p1 [ 8,12 Sar (2,, uxx)) ( (a, (x))))
J8
(33)
Then we repeat this for other derivatives and we end up with the following expression:
(-1)%zt.... ..... (-1)asM [
hx, u(x)))(() ...
...0%r SMr (()()))(x))dd
(34)
As the integrand zero outside of B, this can be rewritten as:
(-1)larp (2, ax (a)()))()))d
(35)
or more compactly, using the functional defined in Definition 3 as:
F(Ep, Φ.
(36)
Therefore Equation 27 can be written as:
Σ F(Ep, u, 6) // lo(o,, u(=)))(x)) da
(37)
p=1
Thus, we proved that Equation 25 is true and only
Σ F(Ep, u Φ) Φ) - - [g(or, u(x))
(38)
for all testing functions
17

---

Under review as a conference paper at ICLR 2023
B.3 EXAMPLES
The examples of VR-PDES can be found in Table 5
Table Examples equations which are ErtationnllReadee
Name
Equation
Linear VR
Damped wave eq. with a source
Ust + pü, A x g(x)
Gauss law
V E = P/e0
Burger's s equation
" + uul Mllz
Navier-Stokes equations
+ (u V Vju V u 1/pVP +g
Korteweg- De Vries equation
Ug + Uzzz Gull = 0
Saannoto-Sivashisky equation
u1 Uxz Urrzz Wtlzz 0
Fisher's s equation
U fills = ru(1 4))
Liouville's 's equation
Ur + Uyy 0u
Porous medium equation
(²((") =0
Sine-Gordon equation
Ust UgA = sin (t)
C D-CIPHER
C.I REWRITE THE INNER OPTIMIZATION AS CONSTRAINED LEAST SQUARES
Let us rewrite the objective in Equation 13
ΣΣ
d=1 F(Spep al)) 26) g(æ; (() (x))
Σ
(39)
First, let observe that
F(BpEp á (d) ).) = 1 hp(ar, (i(00
= =p-
(40)
(d.a) E R. defined
we 2p
(d,a) / (z) (1)) (o))
20
(41)
52
Moreover, we define (d,,) E as
(d(d))) / g(x; (x))
(42)
we can rewrite expression 39 as
D
S
Σ Ap=1
Bp=(de0)
(d.o) (d.o)
(43)
d=1s=1
Now, the sum over p can be written dot product between (d.,) E RR and B E RRP, We can also combine the sums over d and s. We obtain
Σ [d
(g.2 w(dd)))
(44)
which is exactly the same as the objective in Equation 14
C.2 PSEUDOCODE
The pseudocode of D-CIPHER is presented in Algorithm 1
18

---

Under review as a conference paper at ICLR 2023
Algorithm D-CIPHER
Input: Observed fields D (2
Id=1 grid G
Input: Symbolic regression optimization algorithm
Input: Smoothing algorithm S
Input: Testing functions {
Input: Dictionary Q = {2p (1 Ep (ap, Ap:
Step
Output: Target PDE
á ((v(d) Ad E [D]
Step
initialize matrix Z E DXS x R.
zp zfd..)) + fo hp (xr, â (d) (x)) (0x))())))))
procedure Loss(g)
initialize vector w E R Dx5
w(d.3) +- Ja g(r, ((d) (2))
3 < COLLIE Z, w)
Section 6
- 2 w
return L end procedure O(LSs)
Step
initialize vector w E RDx
(d.s) t Sn (a, ald) (z)dx
B COLLLE w)
Section 6
return Ep (pép(u((x) - g( u(x)) =
Step 1: Choose the dictionary
Dataset (Dlo D
0₂2
u²uu
2uu 0, (u?)
d==1
Test functions
Step 2: Estimate the fields
Compute (Eq. 15)
Z Rs I
[i()) d=
Step 3: Optimization
Compute (Eq. 15)
Symbolic Regression
g: R M+N
R
e.g. log(t) ex (1)
Ros
CoLLie
Loss
min ||B wllz
||
Figure 5: This diagram describes how the algorithm works. After the optimization procedure finished, we get the best found closed-form function and use CoLLie to find the best vector B3. The found equation has the form E (u))) g(x,, u(x)) 0
C.3 DIAGRAM
The diagram in Figure depicts how D-CIPHER works.
19

---

Under review as a conference paper at ICLR 2023
C.4 TESTING FUNCTIONS
Testing functions, by definition, need to be sufficiently smooth functions (in C class) with compact support. Moreover, the result proved in (Qian et 2022) suggests that these functions should be subset of Hilbert basis in L.2 space. In particular, that means they should be orthonormal.
We B-Splines [De Boor 1978) as the testing functions in our experiments because we can control their smoothness and the derivatives are easy to compute. We scale and shift them appropriately so that they are orthonormal.
D COLLIE
D.I LAGRANGIAN
The problem that CoLLie is supposed to solve is constrained least-squares optimization defind as:
minimize ||Az - b|1
(45)
subject to ||z||
where E Rx has a full column rank, E T and = E I for some m, n E N.
We consider the Lagrangian L 12 x R R associated with this problem (Boyd et al. 2004)
defined as
L(z, A) 1|Az - b12 >(11=111- -1
(46)
Now let us define 2: - RR" as
$ (A) arg min L(z A) = arg min |Az b112 + >||2111 =ER
(47)
=
The goal of our algorithm is to find X E R such that |a(A")] = 1. Let us define a function
q: R
IR
(()) = ||(())
(48)
The goal can be phrased as finding R such that ((x)) =
Let us note that (0) is just solution to the ordinary least squares (OLS) problem with no constraints and its norm ((0).
D.2 EXTENDING LARS
Case 1. {(0)= 2
If we assume that is continuous then q is also continuous. From the continuity and the fact that
lomx-++oo ((A) and q(0) we infer that there exists À 2 0 such that (()) 1. Moreover, for À 2 0 the problem in Equation 47 is the same as in LASSO (Tibshirani 1996). Therefore we just need to perform LASSO for different À and choose the one that gives the solution with LI norm equal to 1.
To do in practice we use Least Angle Regression (LARS) (Efron al. 2004 a popular algorithm used to minimize the LASSO objective. It generates complete solution paths, i.e., a function R R defined as
c(A) = arg min A bll >||=||1 ZER"
(49)
which is equivalent to & for À > O. An illustration of LARS solution paths be seen in Figure 6 Each line corresponds to a function e; which describes the coefficient for the 7th covariate. The paths are defined from some X where all (A()) 0 to À where c(0) 2(0). In other words, the solution paths cover the whole range of constraints from the strictest, effectively imposing the L1 norm of to be 0, up no constraints, solving the OLS problem.
The solution paths from the LARS algorithm are piecewise linear and the outputs are the values of the coefficients for points (No ... An where the slopes change. We calculate the norm
20

---

Under review as a conference paper at ICLR 2023
at each of these points, ||c(A,)]l1. and find E [n]] such that ||c(A)-1)|ı < <1 ||<(A,)]lz As each c, is a linear function on { and we know both c(Ay_1) and c(Ay), we can effectively search for À E Aj-1 such that |[((A) = 1. The search can be performed by any root-finding algorithm. We use Brent's method (Brent 2013).
Case 2. 0 (()) 1.
This is much more difficult as it corresponds to solving the problem in Equation 47 for < O. The solutions given by the LARS algorithm are too small. In fact, the solution with the biggest norm is c(O)- = =0),, the OLS solution, with norm exactly (()) < 1.
To address this challenge, we propose the following heuristic. We extend the solution paths gener- ated by LARS beyond À = for À < O. We assume that the paths will continue to be piecewise linear and that they will keep the slope they have in the last interval [An 0, An- Let denote this slope as
Aq (0 (1n-1))
(50)
A=-
This is graphically represented in Figure 6 Formally, these extended paths, č: - R are defined as:
(AA) = [c((), (G(0)) + AdG,
X
(51)
x<0
Now, we want to find À < such that [|ε())]l = 1. To achieve this in practice, we first make the following observations.
For any À < 0 we say that (A) is on the right side if (A) S 0 and we say that (A)) on the 0. In other words, being on the wrong side just means that the path is yet
wrong side if (A) to cross the x-axis if we keep decreasing A. We can easily find X' such that for all x < X all (A) the right side (none of the paths will ever cross the x-axis).
0 c0( l [n] (0) >
X min
(52)
AC
If ||(()) 2 we just need to search the interval [X', 0] for such that ||E())] 1.
If ((X\] < 1 then we need to search x < X. However, by definition, for all À < X, all (A)
are on the right side. That means [[ε()]] as a function of is just linear function on the interval (=o0, X"). To see that, let us observe that
|[(A)) = Σ |2(A)) = Σ sign(E(A))E(())
(53)
1=1
Additionally, for À < X' all (A) are the right side, we have sign(E((A) = sign(Act). We can rewrite [ε(A)]] as:
|[(())]| = Σ (igg(A AAc:)
sign (Ac, (O))
Σ sign (c;))
(54)
sign(Ac))n
Σ Ac| x
Therefore the solution can be found using the following equation
X X
(55)
Ac:
Case 3. q(0) = 0. In that case, we just return a precomputed solution to the problem
minimize
A
_
(56)
subject to
| 1
1
=
21

---

Under review as a conference paper at ICLR 2023
which compute by subdividing the problem into 2" guadratic programs and solving each of them
separately using CVXOPT algorithm (Andersen et al. 2013) as described in Section
(A)
(B)
Caeffiic lents
Coefficients
A3
An
As1
Aa
A
A6
Ag
As1
Figure 6: Panel A shows and example solution paths calculated by the LARS algorithm. Panel shows their extended versions defined in Case 2 in D.2 The x-axis is reversed, À decreases as it moves to the right.
D.3 COMPARISON
We perform a comparison between CoLLic and an algorithm based on CVXOPT (Andersen et al, 2013) described in Section CoLLie uses different procedures depending the problem, thus for a fair comparison we generate an equal number of tests falling under Case and Case (two main cases) described D.2 To achieve that, we generate random entry is sampled from the standard normal distribution. We use m = 1000 and n ranging from 2 to We then generate vector & IR" such that each entry is sampled from uniform distribution on [-0.5, 0.5]. Then & is normalized to have L1 norm equal to 1. We sample a number / from uniform distribution on [-1,3] and multiply & by to obtain
matrix A where each
1.2. By this procedure, we are guaranteed that cases ||='|| < 1 and ||=||| 2 1 are equally likely. We then let Az'. The task then to find z with 1.1 norm equal to 1 that minimizes ||Az b 12. As (0) = || I, Case
and Case 2 are equally likely.
We perform 1000 experiments for each n. As losses for optimal solutions can be on widely different scales, report the relative enror between the loss obtained by CoLLie and the loss obtained by the algorithm based on CVXOPT, which seems to always find the optimal solution.
The time is measured single computer with Intel Core i5-6500 CPU (4 cores) and 16GB of RAM.
E EXPERIMENTS
E.
ABLATED D-CIPHER
The ablated version uses the standard MSE loss with estimated derivatives and thus solves the fol- lowing optimization problem:
min min
Σ (p) g(æe (dd)) (x=))
(57)
| d=
where (dd) the observed field, G is the sampling grid, and & (x))))())) requires derivative estima- tion.
The ablated version uses the same symbolic regression algorithm to search over closed-form g and Co for the inner optimization.
22

---

Under review as a conference paper at ICLR 2023
E.2 IMPLEMENTATION
Step 1. For the homogeneous heat equation and Burgers' equation we use dictionary 2 {u, Oru, d,u, O₂ (u²), Hu, 8 ( For Kuramoto-Sivashinsky equation we use dictionary Q {u, Ocu, U,u, O₂ (u²), Hg,, (u(). Lúu, (u(x), ), O2u, 0% (u))). For the damped and forced har-
monic oscillator we use the dictionary {0,, 0. and for the wave and heat equations, we use (0e. 0z, of, 0,0,, ).).
Step 2. Field estimation is performed using the Gaussian Process Regression from the Python li- brary scikit-learn (Pedregosa et al. (01T). The kernel is chosen to be the RBF kernel (Williams & Rasmussen 2006) with an added Wh White kernel account for noise. The observed field imiiially standardized by subtracting the mean and dividing by the standard deviation. Then the Gaussian- ProcessRegressor is fitted the data. The estimated fields are generated by predicting the values of trained Gaussian Process on a full integration grid and then scaling back to their original range (by multiplying by the standard deviation and adding the mean).
Step 3.
The search over the closed-form expression is performed using the symbolic regression library gplearn (Stephens 2022 ). We use a custom fitness function that solves the inner optimiza- tion problem in Equation 13 This inner optimization is performed by CoLLie (Section 6). The
integration is performed using Riemann sums.
Ablated version of D-CIPHER. The derivative estimation is performed by first fitting a Gaussian process (in the same way as in Step 2) and then using the finite difference to estimate the derivative in one of the coordinates for all points in the sampling grid. To obtain higher-order derivatives, Gaussian process is fitted again and the derivative is once again calculated using the finite difference (possibly a different direction than the first time).
E.3 HYPERPARAMETERS
Gaussian process regression. The kernel parameters of the Gaussian Process are automatically adjusted during training. The default bounds of the length scale of the RBF kernel and the noise level of the White kernel are used, i.e., (le 5, 1e5).
GPlearn. We do not perform parameter tuning for the gplearn library and use the same parameters as in D-CODE (Qian et al. (022) except for the parsimony coefficient and the number of generations.
Table 6: Hyperparameters used in gplearn
Hyperparameter
Value
population size
15000
tournament size
20
p crossover
0.6903
p subtree mutation
0.1330
P hoist mutation
0.0361
P point mutation
0.0905
generations
20 and 30
The number of generations is chosen to be 30 for the damped and forced harmonic oscillator and 20 for the inhomogeneous heat and wave equations.
Please check (Stephens 2022) ) for the detailed description of these parameters.
We modify the implementation of the parsimony coefficient. The standard implementation adds to the loss the length of the equation multiplied by the parsimony coefficient. In our implementation, we increase the loss by the parsimony coefficient. This modification is performed because for differ- ent experiments
record the loss widely different scales. To prevent tuning this parameter for every experimental setting we introduce penalty that can work on different scales. The parsimony coefficient is chosen manually by performing experiments for few values. The value used in the experiments is 0.05.
The set of allowed mathematical operations is: (+, -,,,,, x, + sin, exp, log}
23

---

Under review as a conference paper at ICLR 2023
We want to emphasize that we use the same configuration of gplearn in D-CIPHER and its ablated version.
Integration and number of testing functions. For the damped and forced harmonic oscillator we use 10 testing functions and the integration step 0.01. For the inhomogeneous heat and wave equations 100 testing functions and integrate on grid with steps öt 0.01 and ōz = 0.01.
Derivative estimation in the ablated version of D-CIPHER. The Gaussian process is configured the same way described above. The interval used in the finite difference method to estimate the derivative was chosen to be: 10-3
E.4 CHOICE OF EQUATIONS
Equations used in Section 7.1 are canonical equations from physics that often appear in other works about PDE discovery (Rudy et al.] (017). The homogeneous heat equation is a second-order PDE that models how heat diffuses through region. It contains the dissipative term 8.u. Burgers' equation is second-order PDJ used, for instance, in fluid mechanics or nonlinear acoustics (Crighton) 1979). contains the advection term ud,u u and the diffusion term that prevents shock formation. Kuramoto-Sivashinsky equation is a fourth-order PDE used in modelling reaction-diffusiio systems
(uramoto 1980) and is known for its chaotic behavior (Hyman & Nicolaenko 1986).
In Section 7.2 we chose equations of physical significance that have an interesting Ô-free part that is not a linear combination (as discussed in Section ( That makes them impossible to discover by the current methods. A forced and damped harmonic oscillator is second-order ODE. Although D-CODE (Qian et al. 2022) can discover any closed-form first-order ODE, it cannot be used to dis- cover second-order ODEs. Thus there is currently algorithm capable of discovering this equation. Inhomogeneous heat and wave equations are second-order PDEs where the θ-free part is a source (of heat and wave respectively). Moreover, the wave equation does not have the standard evolution form, as it does not involve the å₂ term. Thus even without the source term, most of the current methods cannot be applied directly discover this equation.
E.5 DATA GENERATION
Homogeneous heat equation. The fields were generated by solving the equation 8,u 0,,0g 0 (01 = 0.25) with Neumann boundary conditions 0, (t. 0) (t, X) = 0 and an initial condition u(0, 2) = (z),, where ug is randomly sampled from a Gaussian process. The equation is solved using the implicit BTCS scheme (Kereyu & Gofe 2016) with steps 5t = 0.001 and dz = 0.001. The observed field is generated by sampling (t, x) 10 7 [0, X). evaluating the true field =t, and adding Gaussian noise. T = 2 and X = 2 are used in the experiments.
Burger's equation. The fields with an initial condition u(0, x) = (Ix), where up is randomly sampled from a Gaussian process,
computed by solving d.u + ud,u - 0,0 0 (0, 0.2)
and with Dirichlet boundary conditions uu(t, 0)= = un((). utt,, X) (x) The equation is solved the using Nicnl-Nioosonn scheme (Wani & Thakar 2013) with steps öt = 0.002 and dz = 0.002.
The observed field is generated by sampling (f) E [0. x X], evaluating the field u(t, x) and adding Gaussian noise. T = 2 and X = 2 are used in the experiments.
equation. The solution is the same as the one used in (Rudy et al. 2017 The observed field is generated by sampling (t, .) E [o. TI x [0, X), evaluating the field u(t,, x) and adding Gaussian noise. T = 100 and X = 100 are used the experiments.
Muramoto-Sivashiin
Damped and forced harmonic oscillator. The true fields are created by analytically solving the equation Lf(() + 20,00 (22,v(() + Oz = 0g sin(!4!), where 0, = 0.5, 02 = 4.0, 0g = 5.0, 04 = 3.0, with random initial conditions for u(0) and @,)( The observed fields are then created by sampling E [o, TJ. evaluating ut),, and adding Gaussian noise. T = was used in the experiments.
Inhomogeneous heat equation. The true fields are computed by solving d,u(t, -0 010 x)
fgelat, where 02 = 0.25, 022 = 1.25, 0g = 1.8, with Neumann boundary conditions Ozt(t, 0) = O,u(t, X) 0 and an initial condition u(0, .² up (z),,, where uo is randomly sampled from
Gaussian process. The equation is solved using the implicit BTCS scheme (Kereyu & Gofe 2016)
with steps 5t 0.001 and dæ 0.001. The observed field is generated by sampling (1, x)
24

---

Under review as a conference paper at ICLR 2023
[0, T [0, X], evaluating the true field u(t, x) and adding Gaussian noise. T 2 and X 2 are used in the experiments.
Inhomogeneous wave equation. The true fields are computed by solving -f
sin(t)), where 0 1.0, 02 = 2.0, 0g = 3.0, with Dirichlet boundary conditions u(t, 0) (o()), u(t, X)
(.).. where uo is randomly sampled from Gaussian process and specifies the initial condition u(0, 2: = (x()). The equation is solved using the Implicit Difference Method
(Mitchell & Griffiths 1980) with steps dt 0.001 and dæ 0.001 The observed field was
generated by sampling (1, xt E [0, 7 x [0, X]. evaluating the true field u(t, x) and adding Gaussian
noise. T and X were used in the experiments.
E.6 EXPERIMENTAL SETTINGS
Homogeneous heat equation
Noise ratio (σR): 0.001, 0.01, 0.1 Number of samples (D): 10
Domain (9)) [0, 2 [0.2]
Grid (9): {0, 0.07...... 2} x
Burgers' equation
Noise ratio (ØR): 0.001, 0.01, 0.1
Number of samples (D): 10
Domain (?): [0,2] [0.2]
Grid (9):
2} x {0,0.l 00,......... 2}
Kuramoto-Sivashine equation
Noise ratio (ax): 0.001, 0.01, 0.1
Number of samples (D): 1
Domain (2)): [o, 100] x [0.100]
Grid (G): {0, 1.67, .... 100} {0, 1.67, 100}
Damped and forced harmonic oscillator Default values are in bold.
Noise ratio (ØR): 0.001, 0.005, 0.01, 0.1, 0.2, 0.5
Number of samples (D): 1, 2, 5, 10, 15
Domain (9)) [0,2]
Grid (g): f0, 0.08. 2), 00.. Inhomogeneous heat equation
,2), (0,02.. ...,, ...,, 00.0.4.....
Noise ratio (an):: 0.05, 0.1, 0.2
Number of samples (D): 10
Domain (?): [0.2] [0.2]
Grid (9): {0, 0.07. ...... 2} x (0,0.07, ..,))
Inhomogeneous wave equation
Noise ratio (ag): 0.001, 0.01, 0.015
Number of samples (D): 10
Domain (R): [0,2] [0.2]
25

---

Under review as a conference paper at ICLR 2023
Grid (G): 10,0.07..... ....... x (0,0.07, ..,,,,
E.7 BENCHMARKS
An important hyperparameter in PDE-FIND and WSINDY is the library θ used. We impose that Q and θ {8,u) have the same number of elements. Moreover, solutions given by PDE-FIND and WSINDy are scaled to have the L1 norm equal to 1. Both of these measures are undertaken to ensure that RMSE error is comparable between the algorithms.
For the homogeneous heat equation and Burgers' equation we use a library
O {u, 02 Lu, ud,u, u.
For the Kuramoto- -iinnooo-ivashiskky equation, we use a library
θ {u, 0z Lu, S.u, Hu, ud,u, u, udzu, udz, v02).
In the experiments, we have not optimized for the derivative-free part as it is identically equal to in all equations.
E.8 CORRECT FUNCTIONAL FORM
To measure success probability we need to establish whether two closed-form functions match. The previous approach (Qian al., 2022) considered their functional forms, i.e., expressions where all numeric constants replaced by placeholders. By this measure, functions sin(3x) and (3. match as they have the same functional form sin(Cr), where C is a placeholder.
However, this definition is quite restrictive because functions sin(3x), sin(3x) + 0.001, 1.001 sin(3z), and sin(3z + 0.001) all have different functional forms.
We consider it open challenge to design good metric that would meaningfully reflect whether the correct equation is discovered. We propose the following.
For a target function f. we consider its augmented form J. defined as f(x) = CC C4) + C2. where all C placeholders. Then all numeric constants are turned into placeholders well. In the end, we combine the constants. For instance, C1 C2 becomes just C3.
As example, let consider a function fzz 1x The augmented functional form is created in the following way:
1. Augment: C1 x 2x + C2
2. Replace: C₁ x Cx + C2
3. Combine: C Cgx+Ca
We perform this procedure for the target function. We now take the standard functional form of the candidate function and check whether it matches the augmented functional form of the tar- get function, taking into account that some of the constants might not be present in the candidate expression.
To aid in this procedure, we use Python library for symbolic mathematics. SymPy (Meurer al, 2017).
E.9 COMPUTATION TIME
The average computation time for a single experiment with the damped and forced harmonic oscil- lator is 281 seconds with standard eиror of 4.5 seconds. The average computation time for single experiment with an inhomogeneous heat equation is 68 minutes with a standard error 38 seconds. This time is measured single computer with Intel Core i5-6500 CPU (4 cores) and 16GB of RAM.
The experiments time for all experiments (all seeds, all equations, all experimental settings, and both versions of D-CIPHER) is 65 hours.
simultaneously on computers like the described above. The total
26

---

Under review as a conference paper at ICLR 2023
E.10 LICENSES
The licenses of the software used in this work are presented in Table 7
Table 7: Software used and their licenses
Software
License
gplearn
BSD 3-Clause "New" or "Revised" License
evxopt
GNU General Public License
cvxpy
Apache License
sympy
New BSD License
scikit-leam
BSD 3-Clause "New" or "Revised" License
numpy
liberal BSD license
pandas
BSD 3-Clause "New" or "Revised" License
scipy
liberal BSD license
python
Zero-Clause BSD) license
ADDITIONAL EXPERIMENTS AND DISCUSSION
F.1 SAAAPP-LOTKAA-MCEENDEE MODEL
We test D-CIPHER population model called the Shanpe-Lokaa-Meenddrck model Webb (985). The model deseribed by the following equation
0, a) + dau(t, a) + m(a)u(t,, a) = 0
(58)
where m(a) is age-specific mortality rate. We choose m(a) 2e We 1.5 in the experiments.
We note that the derivative-free part of the target PDE cannot be expressed as a linear combina- tion of functions from finite dictionary if the parameters are not known priori. D-CIPHER is uniquely positioned among other discovery algorithms as the only technique that can recover any closed-form expression. We show comparison between
mortality rate that can be represented D-CIPHER and the Ablated D-CIPHER in Table
model. We report the success proba- bility of discovering the Ò-free part and the Average RMSE of the Ô-bound part. Standard deviations are shown in brackets
Table 8: Simulation results for the
Sucess Probability
Average R.MSE)
Method
OR = 0.001
0.01
OR = 0.001
0.01
D-CIPHER
0.6 (0.15)
0.5 (0.16)
0.007 (0.0008)
0.008 (0.0011)
Abl. D-CIPHER
0.2 (0.13)
0.2 (0.13)
0.017 (0.0009)
0.017 (0.0008)
F.2 DISCOVERING SYSTEMS
Discovering systems of PDEs is much harder problem than discovering single PDE. One of the issues is the fact that we would call indeterminism. It follows from the following fact. If a vector field u is solution two differential equations fi and fz, i.e.,
(x (4) (x), @lk| u (d) (x)) Vo E Ω
(59)
fz u (z). l u ( (z)) V' E Ω
then it is also a solution to any linear combination of these equations, i.c..
A1 S1 (xr u (d) (z:), gl n (4) (a:) (x))) da f2(Ø, u (d) (m)) gl (x)) Va E Ω
(60)
27

---

Under review as a conference paper at ICLR 2023
for any A1, A2 R. Moreover, equations can sometimes be differentiated to yield equations. Let us take as an example the Cauchy- Riemann equations defined as:
i - 42 = 00, + 2 = 0
(61)
Let us assume that we have a true vector field (u1. Uz) that satisfies both equations. following equations also satisfied
Then the
0_ 02 0,40 + 0yu1
(62)
0r 02 0z42 - 0, =
We can also differentiate the Cauchy-Riemann equations to arrive at the Laplacian equations for and u2.
62 0fu1
(63)
8u + 02 = 0
We could also combine the first-order equations with second-order equations consider even higher-order derivatives. Although all these equations are compatible with our vector field (u1, u2), not all of them are equally desirable to discover. That is why we believe that in any algorithm for discovering systems of differential equations substantial expert knowledge or inductive biases have to be encoded to guide the algorithm into the right equations.
Current methods do not consider systems of equations or consider system of equations of a very particular form. In the latter case, each equation models a derivative with respect time of different scalar field. The system is assumed to look like this:
uuu fi u(x), al u(x))
0,ua = f2 uax), glki u(+))
(64)
0ug fL u(x), alkl u(x))
In addition, the LHS of these equations is often assumed only contain spatial derivatives.
D-CIPHER can be used to discover some systems of differential equations if enough prior knowl- edge is provided in the choice of the dictionary Q. Moreover, the discovered equations are not required to have a particular evolution form as is the case in current approaches.
Firstly, we note that D-CODE Qian al. 2022) has been shown discover a system of equations that looks like Equations (4 when all equations are first-order ODEs. As D-CIPHER reduces to D-CODE when applied to first-order ODEs, it also capable of discovering such a system.
We demonstrate that D- D-CIPHER is able to discover both Cauchy-Riemann equations if we use two different dictionaries. Each of the dictionaries yields a different equation. Additionally, well-known fact that if a vector field satisfies Cauchy-Riemann Riemann equations then the constituent scalar fields are harmonic, i.c., they satisfy Laplace's 's equation. Based on the same dataset we are also able to discover both Laplace's equations given another set of two different dictionaries. We note that Laplace's 's equation does not contain O term, so most of the current methods cannot be directly applied to discover this equation. The results are presented in Figure 7 The dictionaries used to discover each of the equations are the following.
For 0,u1 0,u22 0
Qi {0.m, Uy Le, 0,
For 0g 0, 0 we use Q2 {0,u2, 0 11 0 u2 02
For 0.u1 0uu1 0 use Q3 [0.41, 0 u1 0z 0g un}
For 02 u2 0 Q4 (2 U2 02u2, 0f }
In the experiments, we have not optimized for the derivative-oo free part as it is identically equal to in all equations.
28

---

Under review as a conference paper at ICLR 2023
8y41 8,uu 0
8,42 + 8,01
0.06
0.06
Avg. RMSE
0.04
Avg. RMSE
0.04
0.02
0.02
0.0
10
10-4 Noise ratio (a,)
10
0.0
10
10-2 Noise ratio (a))
10
D-CIPHER
8²U1 0²u; 0
8² + 80 =0
0.24
Avg. RMSE
Avg. RMSE
0.16
0.3
0.08
0.0
10
10
10
10
10 Noise ratio (og)
10
Noise ratio
Figure 7: Simulation results for both Cauchy-Riemam Riemann equations and two Laplace's equations. We report the Average RMSE of the Ö-bound part in different noise settings
F.3 INCREASING THE SIZE OF THE DICTIONARY Q
We test D-CIPHER the Sharp--Lokk-MMKenddrin model (Equation 58
with different dictionar-
ies. We start with a small dictionary Q {0,u, d.u). and we create every new dictionary from the previous one by adding one more extended derivative. The final dictionary contains 10 elements, O10 = {a,u, 0,u, Lu, 0fu, uf 0,0,u, u, 0, (u²), (u²). (z(2), (a(²),,f((") 0f (1²",, Og (ur"). The Average RMSE the 0-bound part shown in Figure g We do observe any increase average error. Note that in these experiments we just focus on the Ö-bound part and do not optimize the Ò-free part.
0.015
Avg. RMSE
0.010
D-CIPHER
0.005
2
Size of the dictionary Q
Figure 8: Simulation results for Sharpe-Lotka-MeKendee model. We report the Average RMSE of the 0-bound part for different sizes of the dictionary Q
F.4 COMPUTATIONAL COMPLEXITY
We want to emphasize that PDE discovery is not time-critical application (usually this process is performed manually by scientists) and we believe D-CIPHER's computation time is acceptable for such a task. In this section, we describe which parts of the algorithms are most computationally intensive.
29

---

Under review as a conference paper at ICLR 2023
Computation in D-CIPHER is performed in Step 2 and Step 3.
Step 2. Computational complexity of Step depends on the choice of the smoothing algorithm. We want to emphasize that the user can use any smoothing algorithm based on their domain knowledge and experience, including spline regression, LOWESS, and Kalman filters. Gaussian Process has time complexity (n²) where is the number of data points in grid G. D-CIPHER is specifically designed to work for sparse and noisy data, so we have not encountered major computational issues while performing Gaussian process regression. We also note that significant progress has been made in adapting Gaussian processes for datasets with many data points Liu et al. (2020).
Step 3. D-CIPHER consists of two optimization loops. The outer optimization is performed by symbolic regression algorithm (in our case genetic programming). The inner optimization is per- formed by CoLLie (Section(6). Searching through space of closed-form expression requires testing many candidate equations. D.CIPHER is designed to work with many different algorithms for sym- bolic regression and it is advised to choose an algorithm that can search through this space most efficiently.
CoLLie was specifically designed to solve the optimization problem as quickly as possible with a minor accuracy trade-off (see Figure 1. It is based on LARS which has time complexity O(m)²) Efron et al. (2004) where m is the number of samples and n is the number of features. case, P is usually small as it corresponds to the size of the dictionary and m SD. In our experiments, S was set up to 100 and D up to 10. Overall LARS is performed very quickly. The additional steps in CoLLie require only few arithmetical operations (linear in n) and a possible root searching of a single variable function that is efficiently implemented using Brentq algorithm Brent (0133).
our
Other important parts of the algorithm are the numerical integrations. One such integration is per- formed at the beginning of the algorithm to compute the matrix Z (see Algorithm D. It does not contribute much to the computation time as it is performed only once. The other integration is performed for each candidate equation to compute vector w. Also, substantial time is spent on computing the values of the candidate function g used in the integration. Fortunately, both of these operations can be implemented as vectorized operations which are designed run very efficiently on modern hardware. We also discourage overly long equations for (check the discussion in Ap- pendix E.3) to limit the number of operations performed.
F.5 CHALLENGES OF DERIVATIVE ESTIMATION
One of the advantages of D-CIPHER compared to other methods is the use of the variational for- mulation of PDEs that allows it to circumvent derivative estimation. This is important derivative estimation is challenging, especially in noisy settings with infrequent sampling. The problem be- comes more pronounced the higher the order of the derivative. To demonstrate these issues, perform a series of synthetic experiments.
Qualitative study. First we qualitatively show how challenging the task of derivative estimation is. We generate an observed trajectory for the damped and forced harmonic oscillator. Then we estimate this trajectory using both Guassian Process regression and Spline regression. As shown in Figure g (Panel A), the estimated trajectories are very close to the true trajectory. Then we estimate the first derivative (Panel B) and the second derivative (Panel C). We show the standard finite difference methods as well as derivative estimation techniques using Spline regression and Gaussian Process regression. In both cases we see that the estimated derivatives do not match the ground truth (calculated analytically) as closely as in Panel A. Moreover the mismatch for the second derivative seems to be bigger than for the first derivative.
Quantitative study. We investigate this relation quantitatively for the damped and forced harmonic oscillator and the wave equation. For the oscillator we generate an observed trajectory and then we estimate its derivatives, up to the fourth order using both finite difference and Gaussian Processes. We then compare the derivatives with the analytically calculated ground truths and measure root mean squared error. The results are shown in Figure TO (Panel A). We can see that the error increases the higher the order of the derivative. For the wave equation, we perform a similar experiment but this time we estimate different mixed derivatives. We consider any mixed derivative 0f where
30

---

Under review as a conference paper at ICLR 2023
(A)
(8)
GG)
0.4
0.2
-2.5
Ground truth Firite a fleeecce 5p
Observed dbta
0.0
1.0
4.2
4.4 0.0
-1.5
Figure 9: Panel A shows that estimation of the true trajectory be performed successfully by both Spline regression (Spline) and Gaussian Process regression (GP). Panel shows the estimated first derivative and Panel C shows the estimated second derivative. We observe that the higher the order the less accurate the estimate.
i.j (0,1,2), We demonstrate the results in Figure 10 (Panel B). increases the higher the order the derivative (i + j).
We observe that the error
(A)
(B) a
a
Finite Difference
250
GP
80
0.01 0...0.0.0..0222222 0.49
200
RMSE
150
0.05=0.020.22 1.64=0.73
100
50
a² 0.110221.670000 8.62+2.99 2
0
1
2
3
Derivative order
Figure 10: Panel A demonstrates the eπтог between the estimated derivative and the ground truth increases with the order of the derivative (performed for the damped and forced harmonic oscillator). Panel shows that the same happens for the wave equation. The (i.j) entry of the heatmap should be interpreted the RMSE between the estimated derivative 0,09 and the ground truth.
F.6 CHALLENGES OF PDE DISCOVERY AND HOW WE ADDRESS THEM
PDE discovery a very difficult task with many challenges. Table 9 we summarize some of them and describe how our work addresses them.
Table 9: Some challenges of PDE discovery and how address them
Challenge
Novel contribution (how we address the challenge)
Noisy measurements Large diversity of PDEs
We use the variational loss function (Equation 10
We allow any closed-form 0-free part and require evolu- tion assumption (Section 2
Learning compact equations
We penalize long Ö-free parts and use LI normalization for the d-bound part (Section [5] Appendix [E.3
Efficient search
We develop a quick algorithm, CoLLie, for the inner opti- mization (Section O
Chaotic systems
We do not estimate the initial conditions (Chen al. 2018
perform forward time stepping (Long et 2019 which are computationally unstable for chaotic systems
31

---

Under review as a conference paper at ICLR 2023
F.7 SIGNIFICANCE OF THE NEW NOTIONS
In this section, we want to justify that the new notions we introduce (evolution assumption, lin- ear combination form, derivative-boune part, derivative-fee part, Variational- Ready PDEs) are im- portant theoretical contributions that help us understand the landscape of different PDEs from the machine learning perspective.
The definitions we introduce let us characterize different classes of PDEs. These new notions com- plement the standard recognized PDE classes such as semi- aeni-linear, quasilinear, hyperbolic, etc. These standard classes were introduced predominantly to characterize the solving techniques or the prop- erties of the solutions, whereas the notions we introduce relate to the difficulty of discovering such equations from data.
Significance of linear combination form and the evolution assumption. Table T0 we demon- strate how the presence of the two assumptions, the linear combination form (LC) and the evolution assumption (EA), influences the optimization problem. We see that with both assumptions, the problem is relatively straightforward and reduces sparse linear regression. With only one of the assumptions present the problem becomes more difficult. With neither of these assumptions, the problem becomes very difficult and requires some other assumptions. D-CIPHER does not make either of these assumptions but it assumes the PDE to be of the form described by Equation T2
Significance of O-bound part, Ò-free part and Variational-Ready PDEs. The difficulty of deriva- tive estimation has been of the main challenges of PDE discovery. The variational formulation allows to circumvent derivative estimation and thus more robust noisy data. Previously, the variational formulation has been applied only to a subset of equations in linear combination form and with the evolution assumption. We observed that any restrictions that the variational formu- lation might put on the equation come from the terms containing the derivatives. Thus we define the derivative-buund part and the derivative-fie variational formulation. That allows us to define Variational- Ready PDEs as currently the broadest an important contribution as
part of the PDE due to their significance for the
class of PDEs that admit the variational formulation. We believe methods requiring derivative estimation underperform in settings with high noise. This definition outlines the current limits of any method that circumvents derivative estimation in that way.
Table 10: This table demonstrates how the presence of the two assumptions, the linear combination form (LC) and the evolution assumption (EA), influences the optimization problem. References: [1]
Brunton et al. (2016). [2] Rudy et al.] (2017)((3 Kaheman et al. (2020), [4] Qian et al. (2022)
LC EA Equation form
Optimization problem
Examples
Yes Yes 0.u3- = Ep= Oplp(a, alk) (l
Relatively easy. Can be formulated as finding sparse solution to linear least squares (similar to ridge regression)
[1,2]
Yes No Ep Opfn(x, u(x), al))(()))
Medium difficulty. Can be formulated as P sep- arate linear least squares problems as above Medium difficulty. Find [4] g using symbolic regres- sion
[3]
No Yes dru, = g(x, u(x))
No No (l = -0
Very difficult. Requires other assumptions (in this work, Equation 12
D-CIPHER
F.8 COMPARISON WITH QIAN ET AL. (002
To clarify our contributions, we compare D-CIPHER to Qian et al. (2022) which we think is closest in spirit to the algorithm we developed both of them allow any closed-form derivative-free part
32

---

Under review as a conference paper at ICLR 2023
and variational formulation circumvent derivative estimation. We also note how the new notions we introduce help us compare the two works.
Algorithm presented in Qian al. (2022) discover any first-order explicit closed-form ODE, ic., an equation of the form.
(tu = g(
where g is a closed-form function. It uses the variational formulation ODEs to circumvent deriva- tive estimation.
Table II: Comparison between D-CIPHER and Qian et al. (0022)
Property
Qian et al. 2022)
D-CIPHER
Applicable PDEs
No
Yes
Higher-order derivatives
No (only first-order)
Yes
Requires evolution assumption
Yes
No
Derivative-bound part
Fixed: 0,uy
Learnd: Epe 3pE
Derivative-free part
Any closed-form function No
Any closed-form function No
Derivative estimation
Qian (2022) can be considered special case of D-CIPHER. We can recover it from D-CIPHER by choosing the dictionary to contain only one element, i.c., 2 = drug.
As the derivative-bound part is fixed, every ODE of that form admits the variational formulation. This is not true for PDEs as there are derivative- bound parts that might prohibit the variational formulation. Thus D-CIPHER required careful consideration the appropriate class of equations to search over.
As D-CIPHER needs find both the O-free part (function g) and the 0-bound part, the optimiza- tion problem is much more complicated. That is why we restrict the derivative-bound part of the PDE be spanned by terms from the pre-specified dictionary and develop an efficient optimization algorithm, CoLLie. We emphasize that, as is the case for Qian et al. (0022), we do not put any constraints on the derivative-free part the PDE, apart from being closed-form.
F.9 ERROR BOUNDS
While we would like to have error bounds for the discovered PDE, we note that the problem we solve is significantly more difficult than the one considered in other works. The space of PDEs we consider is much more complex than the space of PDEs in a linear combination form. other works (e.g.. Rudy et al. (0177), Messenger & Bortz a the PDE is basically vector in Rº and the discovery task is mostly reduced to finding sparse enough vector that approximately solves certain linear equation. Of course, there is lot of literature that aids in establishing error bounds in such problem settings, However, D-CIPHER searches over a space Rº x CF CFE + +N), where CFFE(++H NN) V) is space of closed-form expressions in M+ variables. This space is combinatorial in the functional form and continuous in real constants. This makes it very challenging to derive any error bounds.
Hamogeneous heat equation
Burgers' equation
Aaramoto-SSvashiig equation
0.25
Reg RMSE
0.2
Ivo PMMSIE
GP
0.15
0.15
Nearest
0.1
0.2
Linear
0.05
0.11
0.05
Cubic
0.0:
10
10-2
10
16-2
10
10
10
16
Noise ratio (as)
Noise ratio (ax)
Noise ratio (mx)
Figure 11: Comparison of different estimation algorithms that can be used in D-CIPHER. GP Gaussian Process regression, Nearest Nearest point interpolation, Linear Linear interpolation, Cubic Cubic interpolation.
33

---

Under review as a conference paper at ICLR 2023
F.10 COMPARISON BETWEEN DIFFERENT ESTIMATION ALGORITHMS
Step 2 of D-CIPHER requires estimating the fields. We emphasize that any choice of reconstruction algorithm can be used, and it should be chosen based the application and domain knowledge. In our experiments, D-CIPHER is implemented using Gaussian Process regression (Williams & Ras- mussen 2006). In this section, we investigate other common interpolation algorithms. We imple- ment D- D-CIHHER with different estimation algorithms. In particular, we compare Gaussian Process (GP) against: Nearest point interpolation (Nearest), Linear interpolation, and Cubic interpolation (Cubic). The implementation details of these three algorithms can be found in the scipy (Virtanen) (ct al. 2020 documentation. The results are presented in Figure [T algorithms that produce smoother functions (GP, Cubic) tend to give better results.
see that the estimation
34


Reviews:
[{"review_id": "g0fIK3veCF", "reviewer": "Reviewer_HGvL", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "Authors propose an algorithm called  D-CIPHER to elicit differential  equations from data.  Authors also propose a new optimization procedure called CoLLie, to help D-CIPHER. Empirical examples are provided. \n", "summary_of_the_review": "This is an empirical study on an algorithm that is claimed to elicit PDEs. However, the lack of a proper theoretical framework and results and a lack of discussion about caveats and limitations make the work weak.\n\n\nAfter author comments, several issues are clarified. While some doubts remain, I am changing my scores to better than earlier positions.", "strength_and_weaknesses": "Strengths:\n\n\nAuthors have a good grasp of some of the challenges of eliciting differential equations from data. They have also planned the scope of the problem carefully, as can be evidenced from the notions like such as evolution form, evolution assumption, derivative-bound part, and derivative-free part, variational-ready PDE. \n\nWeakness:\n\nThe paper is one of bunch of definitions and some empirical analysis. It is not clear that the algorithm D-CIPHER can get past some of the major challenges with eliciting PDEs from data. Essentially, PDE estimation depends on the smooth function used to interpolate between observed data, and hence effectively depends on uncountable imputations between observations. Consequently, the estimated PDE (by any algorithm) tends to be heavily assumption dependent and extremely sensitive to noise. This issue has been deeply studied in the applied math literature and related machine learning literature, but the authors do not seem to be aware of this. The details in the supplementary materials also do not provide any theoretical results or hint that the algorithm indeed provides a consistent estimate of a PDE. I am surprised to find no discussion on or constraints and caveats based on signal to noise ratio, or chaotic systems and other well-known difficult/impossible cases of PDEs. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "eaNeRUoK2h", "reviewer": "Reviewer_54d7", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "The authors propose a symbolic-regression encapsulating approach to discovery of partial differential equations and a heuristic approximate solver of the constrained least-squares problem arising from the definition of learning PDEs as an optimization problem. The approach is evaluated and compared with related approaches on 7 synthetic experiments.", "summary_of_the_review": "The work addresses an important problem in an area of symbolic regression that requires more attention, as approaches to discovery of PDEs are scarce and mainly developed as extensions to either standard SR methods or methods for learning equations for dynamical systems in the form of ordinary differential equations. I am not an expert in the area of PDEs, but the presented results look convincing. The approach is presented in a manner of framework where different modules can be chosen based on the problem at hand. That being said the effect of the choice of these modules at least in the context of the experimental evaluation in the manuscript deserve a more extensive treatment. ", "strength_and_weaknesses": "Strengths:\n- Well written\n- Complete description of the problem and solution\n- Comparative and in depth evaluation on multiple problems\n\nWeaknesses:\n- High-level abstract description\n- Dependence on field estimation and SR algorithm not discussed sufficiently", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "0fy4_NCuR0O", "reviewer": "Reviewer_HS33", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "This paper proposed a PDE/ODE discovery method that includes: (1) a new designed loss function which depends on the variational form of PDE; and (2) an optimization procedure to efficiently solve a constrained least-squares problem. It can discover PDE with nonlinear partial-free part and high-order ODE. ", "summary_of_the_review": "I have provided detailed evaluation in \"Weakness\". ", "strength_and_weaknesses": "Strength: (1) The designed method has the ability to discover the equation that cannot be represented as a linear combination of derivatives, and the empirical results are promising. (2) Using variational form as the objective in PDE discovery is novel.\n\nWeakness: (1) The relation between the challenges in PDE discovery and the proposed method is not clearly stated. For example, to discover the last four equations in Table 4, why original loss does not work? Is the proposed variational objective function necessary to discover this kind of PDE or just for optimization robustness?  Another example is that it claims the proposed method can find high order differential equations in a robust way. Which techniques can reflect robustness and what makes baseline method unstable? I cannot get why the proposed method has special benefit on learning high-order differential equations neither, because it still depends on the pre-given elements in the dictionary. I feel that there are many claimed points which makes me hard to catch which one is the main point that really works. (2) The presentation is poor. For example, the optimization procedure (CoLLie) is listed as a contribution in the introduction, however, there's no details about the algorithm in section 6. Therefore, I cannot justify whether the proposed procedure is new or not. The figures in experiments are not clearly introduced. For example, what is \"Abi. D-CIPHER\"? Why does some lines in Figure 2 disappear? ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject