D-CIPHER: Discovery of Closed-form Partial Differential Equations
Authors: Krzysztof Kacprzyk
Zhaozhi Qian
Mihaela van der Schaar
Venue: NeurIPS 2023 poster (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 phenomena. Moreover, many of them resolve the equation-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 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.

D-CIPHER: Discovery of Closed-form Partial Differential Equations
Krzysztof Kacprzyk University of Cambridge kk751 177mm.cm.... ac
Zhaozhi Qian University of Cambridge zq2240cam. ac. uk
Mihaela van der Schaar University of Cambridge, The Alan Turing Institute mv472@cam ac uk
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 model and better understand natural phenomena. Discovering these equations directly from data is challenging because it requires modeling relationships between various derivatives that not observed in the data (equation- data mismatch) and 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 phenomena. Moreover, many of them resolve the equation-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 measurement artifacts and uncover new and very general class of differential equations. We further design novel optimization procedure, CoLLie, 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.
Introduction
Scientists have been using mathematical equations describe the world for centuries. particular, closed-form differential equations turned out to be one of the best tools to model physical phenomena. differential equation describes relationship between quantity and derivatives (rates of change); it is called closed-form if this relationship is described by mathematical expression consisting of finite number variables, constants, arithmetic operations, and some well-known functions (e.g., exponent, logarithm, trigonometric functions]] Closed-form differential equations provide general description of reality concise representation that amenable to closer inspection 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 discovering differential equations from data
can only identify Ordinary Differential Equations (ODEs) which evolve only with respect to one variable (usually time). In contrast, many natural phenomena described by equations involving many variables (c.g.. spatial coordinates) called Partial Differential Equations (PDEs). Many equations also involve higher-order derivatives.
Partial and higher-order derivatives. Many algorithms 10 42
Detailed discussion in Appendix A.2
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Derivatives not observed. Discovering differential equations from data is challenging because the derivatives are usually not observed in the dataset (equation-data mismatch (42)) This makes verifying candidate equation a non-trivial task. Most of the methods try to resolve this issue by
estimating the derivatives [T0 47. However, derivative estimation is difficult, especially when the data is sampled infrequently or with high noise [4 35 For an illustrative study, see Appendix F 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 make the evolution assumption (defined and explained later) and assume that the equation can be represented linear combination of prespecified functions and differential operators [10 35]. However, many well-known equations, such as a forced harmonic oscillator or an inhomogeneous wave equation, cannot be represented in that way.
Currently, a few algorithms tackle only some of these challenges. In particular, Weak SINDY able to discover PDEs without estimating the derivative by utilizing variational approach. However, the form of the equation is constrained to be in a form amenable for a sparse regression algorithm. D- CODE [42,, on the other hand, uses variational 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 PDEs. The algorithms that do not require the evolution assumption appeared in [3] and (23 but they require derivative estimation and only consider equations represented as linear combinations of prespecified functions.
35
Contributions. In this work, we develop the Discovery of Closed-form Partial and Higher-order Differential Equations in a Robust Way framework (D-CIPHER) that does not estimate the derivatives, 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-bunnd part, and derivative-free part. We use them 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 (Variational Ready 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 RRadyy PDEs. (Section 6)
In addition to the main contributions above, we also develop a new optimization procedure (CoLLie) to help D-CIPHER search through this space efficiently. (Section 7
2 Preliminaries
In this section, we provide background information about PDEs and their variational formulation.
Notation and definitions. We denote the set {1, 21.....,77 n) as [n] and the set of non-negative integers as No- Throughout this paper we let M, N,, E be some natural numbers and let C R-M be an open set inside R-M. A comprehensive table with all symbols used in this work can be found in Appendix A 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 only one independent variable, usually time. Most methods assume that the ODE is explicit and can be represented as a system of first-order ODEs: (I)
úe (t) fy (t, u(t))
where ùy represents the derivative of uj. Then the discovery problem reduced deciding the order of the derivative (usually first or second) and the discovery of f,-
For PDEs, it is not enough to talk about the derivative, as can take derivatives with respect to different variables. We denote the mixed derivative as 80 where a E Na is called a multi-index, and define it as 0 = 0f 0g2 J.. 0i Each 8.:: Oo /ax. tt is a th -oddr partial derivative with respect to D; (the th independent variable
We define the order of a as la = M ag We call dm
non-trivial if [x| > 0.
Throughout this work we assume that the functions we use are smooth enough for the equality of mixed partials (50) to hold. In that case, any mixed derivative can be uniquely specified by indii--ieexx

---

A PDE of order K' is any equation of the form
f (x)), k u (x)) = Va E Ω
(2)
where @lK are all non-trivial mixed derivatives all uy (i [N) up to the K order. call PDE closed-form if J is closed-form.
Variational formulation (VF) of PDEs is way describe PDEs without referring their derivatives [9]. It works as follows: we take a differential equation, we multiply it by special testing function, and integrate. Finally, we perform integration by parts move the derivatives from the dependent variable u onto the testing functions. E.g., for homogeneous heat equation,
0,u -8 + 1 (8.u - u0 dtdr 0 Vo + J/8 - oud dtdr Vi
(3)
For more details, see Appendix A and B By not depending on the derivatives, methods that utilize VF are more robust noise than their derivative-estimatnng counterparts [42 36].
3 Relaxing assumptions while staying robust to noise
Below, we introduce the evolution assumption (EA) and the linear combination (LC) assumption made by the current discovery methods.
Evolution Assumption. Although there is no generally accepted notion of an explicit PDE (as is the case for ODEs), we define evolution form of PDE be an equation of the form
8 (x) f(z, u(x), a [K///a u(x) Væ Ω
(4)
where ||K]/a is []K] with 8 omitted, a is a known multi-index and j E [N]. Note that if M = 1 and |ox| = K then Equation becomes exactly the definition of an explicit ODE,
In fact, many algorithms for PDE discovery assume a particular evolution form 35
We call it an
evolution assumption (EA). However, this assumption requires the knowledge of a and which might not be trivial. Usually, 0 is assumed to be the first derivative with respect to time (8,) [47 but it is not the case for many well-known PDEs such the wave equation or Gauss's law.
Linear combinations. Current PDE discovery algorithms 47 35 12
consider PDEs that are linear in parameters. That means the PDE can be represented as a linear combination of functions, i.e.,
Oppple u(x), a (x))) = 0 V E Ω
where 0, E I for [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 particular, for an unknown 8 E IR functions such as sin (fx:)), ez, or 2:49 cannot be learned by these algorithms.
How to relax these assumptions and still allow for variational formulation? Current methods that utilize VF either assume that the PDE is in an LC form they only work for explicit first-order ODEs. Moreover, all of them also make the evolution assumption. Relaxing L.C is not trivial because not all PDEs admit VF. As in Equation the PDE has be by parts can be performed. Our crucial observation is that for any term that does not contain any derivatives (and thus does not need to be integrated by parts) additional constraints need be put in place. Due to the significance of these terms, we propose the following characterization of PDE.
of terms for which the integration
Derivative-bound part and derivative-fres part. Any PDE can be expressed in the form
fæx u(x), glki u(x)) - g(x, u(x)) = 0 Va E Ω
(6)
where we collect all the terms with the derivatives into f(x, u(x), (z))) and all terms without the derivatives into g(x, u(x)). We call the derivative- bound part and g the derivative-free part (denoted also O-bound and Ò-free). Ò-free part can be evaluated directly given u, whereas the 0-bound part requires access the derivatives. Note that for first-order ODEs, is trivial and equal to út-
Constraints on the derivative-bound part. Although VF does not require any constraints on the O-free part, we still need to put some constraints the d-bound part for the integration by parts to
3

---

Table 1: Columns correspond to challenges outlined in Section I and answer the following questions: Can discover PDEs? Does avoid derivative estimation? Is the evolution assumption unnecessary (Equation ! ? Can it discover any closed-form Ò-free part (Equation ?
Method
PDEs
No & estimation
No evolution assumption
Any closed-form d-free part
SINDY 10
X
X
SINDy-implicit
53
X
X
PDE-FIND 47
/
X
x
x
PDE-Nel 2.0
29
/
X
X
X
WSINDY R5
46
/
V
X
X
D-CODE
4
X
x
D-CIPHER
work. This is what we do in Section 5 where we aim to define currently the broadest form of PDEs that admit the variational formulation using the above characterization.
Optimization challenge. D-CIPHER does not need the evolution assumption and it can even discover some PDEs that cannot be put into the evolution form. Moreover, unlike previous methods, D-CIPHER is not limited to PDEs that can be represented as linear combination of functions (we describe the exact form we assume in Section O This makes the optimization problem much harder as we search over all closed-form functions 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 [2) or they constrain equally both the 0-bound part and Ü-free part to be linear combination of some pre-specified
functions (35) (for more details, see Table I0 and [a Appendix F. One way we address this challenge is by developing a new optimization procedure (Section)
4 Related works
Symbolic Regression. The goal of symbolic regression is to find a closed-form expression that best models the given dataset both terms of accuracy and simplicity. In contrast with the conventional regression analysis which optimizes the parameters of a pre-specified model, symbolic regression aims discovering both the general structure and the parameters of the model. Most of the existing work focuses on developing optimization algorithms. Genetic Programming (26 has been widely
used for that task (49. A different strategy has been employed in AI Feynman 54 55 that uses
neural networks to reduce the search space by identifying simplifying properties like symmetry or separability. Optimization methods based on pre-trained neural networks 20 reinforcement
learning 440] and Meijer-G functions [I 13 have also been proposed.
Data-driven discovery of closed-form differential equations. Data-driven discovery physical laws is an established area of machine learning [6] 44). The pionecring work in that area was SINDy [10 that constrained the space of equations linear combinations of functions from a predefined library and used sparse regression to discover explicit ODEs. It was later extended to include and PDEs 47 Various other extensions were proposed by improving the
implicit ODEs 33 23
derivative estimation and the training procedure [45] 61), adding additional selection criteria [32]
different approach is taken by 29 (an extension of (30)) which uses convolutional and symbolic neural networks. It is important to note that all of these methods still assume the PDE to be a linear combination as discussed in Section 3 (Equation 5) which significantly limits their search space. Some other developments are but they require the exact form of the PDE and only optimize
and learning the library using genetic programming [4 12 62
based on Gaussian processes 44 43
the parameters.
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 [35] 36 46 or closed-form first-order ODEs 42). 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 new and general class of PDEs in Section 5 that admit the variational formulation.

---

Variational-Read, PDEs
In this section, we propose and very general class of PDEEs, the Variational-Ready PDES (VR- PDEs), which can be characterized without referring to the derivative. The VR-PDEs allow arbitrary O-free part but make some minor restrictions the d-bound part. These restrictions allow to use the variational formulation of PDEs to 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 class of PDEs, need the following definition.
Definition (Extended derivative and differential operator). Let a E fal S. K, be a multi-index. Let h BM+N - IR and RM + IR be smooth functions. An extended derivative &, denoted (ox, a, h), maps a vector field u R-M + R to a function Elu] IR - I defined as:
[(u[(z)
( u(x))]
E is called closed-form if and are closed-form. We call & non-degenerate if ||x] 0.
Now, let pE(P) be finite sequence of non-degenerate extended derivatives. The extended differential operator, denoted as Eip] is an operator defined as:
{(p]/u[(x)> = Σ (x) p=1
(8)
Remark. Any linear operator I, = EoEA 0a acting on "j is an extended differential operator. Definition 2 (Variational- Ready PDE). Let Elp] be an extended differential operator, and let g RM+N + R be a continuous function. We denote a Intiaationl-Reada PDE (VR-PDE) by pair (E))) and define it as:
[p]|u](()) g(x, u(x)) O Va E Ω
We extend the standard variational formulation of PDEs (Proposition]] iin Appendix B) from linear PDEs to all VR-PDEs. The following definition is useful in further discussion.
Definition 3. Consider field Ω R&, and an extended derivative & = (a,, a, h). Let
6: Ω + be a testing function (CK function with compact support). We define the functional
u. 0) [ h(w, u (==))((1))))))) (a(z))(<))]dr (())))))
We can now use this functional formulate variational characterization of VR-PDES. Theorem 1. u Ω R,,, where u, E CCK,, is solution to VR-PDE in Equation if and only if
P
F(Ep, u. - [9(2, ((z)) Σ
(10)
for all testing functions @ 0 R.
Proof. Appendix B
This theorem motivates the variational loss function as we expect the left-hand side of Equation 10 to be closer to the closer the canditate PDE is to the true one. To calculate how well a set of vector fields D = (u(d) Id=1 matches VR-PDE (E1P).9) we propose the following loss function.
D 9 ΣΣ d=
L (8
F(E, ul (2) g(22, ul (x)dxr)
(11)
where (0,)5- is set of predefined testing functions.
This novel loss function makes it possible to evaluate to what extent any VR-PDE matches the observed data. This loss can be used as optimization objective in any algorithm that searches over some subspace of closed-form VR-PDEs. We propose D-CIPHER in Section a an example of such an algorithm.
We say Ie
J. is in CK if 8% exists and is continuous for all |a S K.
5

---

6 D-CIPHER
In this section, we formulate the problem of PDE discovery and then introduce novel algorithm (D-CIPHER) to solve it. The diagram and pseudocode are presented in Figure Land in Appendix C
Problem formulation We are given a dataset of fiels D {(u(d) ) d==1 with finite sampling grid 9 C f2 Each o(d) (a:) a noisy measurement, i.e., w G -+ R.V defined as
00 (z) - uj ( + + (x) Vo E G Vjj E [N]
(12)
where (100 (x) is a realization of a zero-mean random variable (noise), each // (/: :2 - IR is CK function, and every true field (d is governed by the same closed-form PDE f. The task is to infer the closed-form PDE, f. from the dataset D {(r(1) 12_1 and the sampling grid G. We assume that f is inside the class of closed-form VR-PDEs (Section 5) and its 0-bound part is inside a subspace of extended differential operators spanned by user- sper-spe dictionary (see Step 1 below).
We propose an algorithm that consists of three steps. the first step, we define the subspace of closed-form VR-PDEs we want to 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
Step 1:
Choose the form and incorporate prior knowledge. A human expert should encode their prior knowledge of the problem into dictionary of -en-eegeeee extended derivatives Q (Ep) We use this dictionary to search over mnite-iimensionll subspace of closed-fore operators spanned by this set. In other words, we assume that the VR-PDE of the form:
P
( - g( u(x)) Va E Ω p=1
(13)
where B E RR,, 99 is any closed-form function of M + N variables, and E,, (Ctp. Ups hp).
For instance, a dictionary might include only the partial derivatives up certain order. For a 1+1 second-order equation that means Q = { Oz atz, 02 } 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 ut, (as in the Burgers' equation), the term O₂ (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 DD consists of noisy and infrequently sampled fields, we first need to estimate the "true" fields (d) from v(d), Any choice of reconstruction algorithm can be used and the user should choose it according to the problem setting and their domain knowledge.
Step Optimize. We minímize the loss function in Equation 11 for the estimated fields (ū ) among all PDEs of the form in Equation] 13 We solve the following optimization problem:
S
1 g(æ, a (x))).
min min ΣΣ d=l
F(Bppp .0)) (0,) -
(14)
As we want to discover both g and (6 we cannot use the standard penalties on B such as the A/ A||||||2 or A||B|1, as the loss would be minimized by O and B = 0. Therefore we put the constraint ||@||] = 1. We choose the norm to encourage sparsity in the coordinates of the vector G..
The inner minimization in Equation 14 can be rewritten as a constrained least-squares problem.
(3 (d.a) w(d.))) 2
min ||| =1
Σ
(15)
[d)
where &,- (App ap, hp) and ((d.)) E R., wd:ss) E IR are defined as
2dd,0)
h (1)) g (00
(16)
(d.s)
g(2, ((() (2))

---

We show the full derivation in Appendix C z (d,s) can be precomputed at the beginning of the algorithm without estimating the derivatives of the reconstructed fields. They can be easily calculated if the derivatives of the testing functions 0, and the derivatives of Op can be analytically computed.
As the optimization problem in Equation 15 has to be solved many times for different closed-form expressions it poses some unique challenges. As standard approaches are not sufficiently fast, we design new heuristic algorithm to solve this problem, CoLLie, and describe the next section.
Step 1: Choose the dictionary
(Ep) p=1
d,u
Dataset (Dla)
ofu
D
0,0,u
dd=1
u 02 (x())
Test functions
Step 2: Estimate the fields
(
Compute (Eq. 16)
Z R R
((((0) D
d=
Step 3: Optimization
Compute (Eq. 16)
g: RR M+N
Symbolic
w RO
Regression
e.g. log(t) ex (i)
CoLLie
Loss
min I|Pl =1
|z
wll2
Figure 1:
This diagram describes how the algorithm works. After the optimization procedure is finished, we get the best found closed-form function and use CoLLie to find the best vector (3. The found equation has the form Ep=1 (p)((()( g(x, uxx)) = 0
7 CoLLie
The problem in Equation [5 from the previous section can be formulated as follows. Given matrix A E _uxn and vector b Rm find vector E R= that minimizes ||A bi such that ||>||| = 1. The task is challenging as the unit L1 sphere is not convex. A method that guarantees an optimal solution is based
observation that the (n 1) dimmeensional L1 sphere consists of 2" (n !)-simplices (which are convex). Minimizing || - bll2 on a simplex is a quadratic program 8 with many available solvers 51 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 heuristic algorithm CoI.Lie (Constrained LI Norm Least Squares) that finds an approximate solution but is significantly faster (Figure 2). We observe that this optimization problem is related to the one encountered in LASSO. Denote zo the solution that minimizes A b) 1/2 (no constraints). If | 2 1, the problem is equivalent to finding À (in the Lagrangian form of LASSO, Equation 47) such that the LASSO solution has the norm 1. Least Angle Regression (LARS) [17 is popular algorithm used to minimize the LASSO objective that computes complete solution paths. These paths show how the coefficients of the solution change as À moves from ( to Xmaz (from no constraints to effectively imposing the solution to be 0). See Figure 6 in Appendix D.2 CoLLie uses these solution paths to calculate the exact solution to optimization problem. The case 0 < ||= < is harder as it corresponds to À < 0. CoLLie addresses this challenge by extending the solution paths generated by LARS beyond
for O. We assume that the paths continue to be piccewise linear and keep their slope (Figure 6 in Appendix
7

---

D.2. CoLLie then uses these assumptions to efficiently find an approximate solution. We provide detailed deseription of CoLLie in Appendix
(C)
Avg. relative error
10
0.009
Relative etror
10
10-1
CoLLie
10-2
0.003
10
CVXOPT
10-3
10 -15
0.0
Figure 2: We compare CoLLie with an algorithm that uses CVXOPT [2) to solve each of the convex subproblems. We report the relative error between the loss obtained by CoLLie and the minimum loss achieved by CVXOPT. Panels B and C show the averages and the distributions of relative errors. The average relative error is below 0.005 and the bulk of the distribution is below 10-7, At the same time CoLLie is orders of magnitude faster (Panel A).
Experiments
We perform a series of synthetic experiments to show how well D-CIPHER is able to discover some well-known differential than current methods when discovering PDEs in linear combination form (Section 8.T). Then we demonstrate it can discover PDEs with closed-form Ò-free part that cannol be expressed as linear combination and thus are beyond the capabilities of current methods (Section (.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.I). For additional information about the experiments (c.g., implementation details, data generation, experimental settings) see Appendix E
First, we demonstrate that D-CIPHER performs better
Table 2: Equations used in the experiments. "L.C" column specifies if the equation can be represented as a linear combination (Equation 5). "VR" column specifies the PDE is Variational-Reaat
Name
Equation
LC
VR
Homogeneous heat equation
8,u 0, 0.
/
/
Burger's equation
8,u + u0zu -0 0103,u =0
/
/
Kuramoto-Sivashinsky equation
0,u + 82 -O+u + uuzu-
V
V
Forced and damped harmonic oscillator
0f + 20, 20,020,u u + 0z O3 sin(Ost)
X
V
SLM model (Appendix F.1
d,u + 0z -F
X
Inhomogeneous heat equation
0,u 01 .& = B2e
X
Inhomogeneous wave equation
0fu 0200u Oze' sin(Ogt)
X
Evaluation metrics. To establish how well a discovered PDE matches the ground truth, we evaluate its ü-free and Ö-bound parts separately. For the Ò-free part, we assign a binary variable indicating whether the correct functional form of the equation was recovered (please check Appendix EE8 for details). For the O-bound part, we measure the RMSE between the found coefficients of / 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 15 as the testing functions and we estimate the fields in Step of D-CIPHER with a Gaussian Process [60], The outer optimization in Step is performed using
modified genetic programming algorithm (26) and the inner optimization by CoLLie (Section 7 We
also show additional experiments with different estimation algorithms Appendix F
8.1 Discovering Linear Combinations: comparison with other methods
We compare D-CIPHER against two variants of PDE-FIND 47 and WSINDY [66 (as implemented
PySINDy library 24 16) with optimization performed by Stepwise Sparse Regression or
All experiment code be found htttttss///gthhhi httpaa//oe./eettttttaaarrr//////

---

Forward Regression Orthogonal Least-Squares [3] We note that D-CIPHER is specifically designed discover PDEs that are beyond the capabilities of current methods, i.e., where the derivative-free part can be any closed-form expression. Current methods are usually tested on equations where the derivative-free part is trivial (identically equal to O). Even though these algorithms are specialized to 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 homogeneous heat equation, and Kuramoto-SSvvshinnky Sivashinsky equation (Figure 3 This demonstrates gain from both the variational loss and the new optimization routine.
Homogeneous heat
Burgers' equation
Kuramoto-Sivashinkky equation
0.5
0.25
Neg MMSE
0.08
0.1
WSiNDy (S)81F)
Ngg RMSE
0.4
Avg. RMSE
0.2
[SBBP]
PDE-FIND (5)
D.CIPHER
0.3
0.15
PDE-FIND
0.04
0.2
0.1
(SW6))
WSINDY (5)
0.02
POE ANN (SIBIF)
0.1
POE FIND (S18IF)
0.05
WSINDY
10
10
0.0
10
10-1
0.0
10
10
Noise ratio
Noise ratio
Noise ratio
Figure 3: Simulation results for the Burgers' equation, homogeneous heat equation, and Ku-
Aonoto-Sivashinsyy equation. We report the Average RMSE of the Ô-bound part of the equation. Note that some of the benchmarks overlap
8.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 [42]. D-CIPHER discovers the correct functional form of the Ò-free part and achieves a low RMSE for the coefficients of f in most of the experimental settings. The performance higher than comparable the ablated version of D-CIPHER, thus demonstrating gain from using the variational approach. We present the results in Figure
Forced damped harmonic oscillator
Probabiity
0.8
Succes Probabiity
0.8
Success Irdaablity
0.8
9.6
0.6
0.4
0.4
0.2
0.2
0.0 30
10-3
10c4
0.0
0.1
0.2
0.1
0.4
0.0
10
15
D-CIPHER
Noiseratio lou)
Sampling interval 18/1
Number of samples (2)
ADI D CIPHER
Rug, RM-SE
0.06
0.12
800
0.04
0.00
Aup.
0.08
0.00
0.02
0.94
10
16
16-4
10*
0.0
0.1
0.2
0.3
0.4
0.0
5
10
15
Noiseratoo (og)
Sampling Interval 1641
Number of samples 1D)
Figure 4: Success probability of discovering the correct θ-free part of the equation and the average RMSE between the recovered d-bound part and the target one across different experimental settings. We compare D-CIPHER against its ablated version (Abl. D-CIPHER).
Inhomogeneous heat equation. D-CIPHER is able to discover the correct equation even in settings with very high noise. performs better than the ablated version, thus showing the importance the variational objective. The result are presented in Table
Table 3: We report the success probability of discovering the P-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
OR = 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 it does 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 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 not far from the ground truth. The ablated version fails to discover the correct functional form and the found 0-free part does not reproduce the correct behavior of the equation.
10
0, 10
1==15 10-
xelsin/2.98c)
x (9))n(2
Error
D.CIPHER
sinte!
2,99 (sinte)) log
Abl. D-CIPHER
Figure 5: 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 x e! sin(3t).
9 Discussion
Applications. As D-CIPHER potentially discover any closed-form Ò-free part, it is especially useful when this part of the PDE captures an essential component of the phenomenon. We demonstrate it by finding the heat and vibration sources as well as the driving force of oscillator. Beyond the spatio-temporal physical equations, D-CIPHER might prove useful in discovering population models structured by age, size, and spatial position (59 58 age-dependent epidemiological models [21],
and predator-prey models with age-structure 41 All these equations are VR-PDEs where the Ò-free parts crucial elements of the equations signifying the rates of mortality, infection, recovery, or growth. We believe that discovering closed-form equations for these systems would prove invaluable understanding their behavior.
Limitations and open challenges. D-CIPHER may fail in some scenarios, either due to challenging experimental settings 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 PDE outside of the VR-PDE class or a O-free part with a complex expression that is difficult to find. We note that we address some of these challenges by utilizing variational approach, defining VR-PDES to be very general class of equations, and designing CoI.Lie, enabling thorough search across closed-form expressions.
Ethics Statement. We want to emphasize that D-CIPHER designed to facilitate the process of scientific discovery by extracting closed-form PDEs from data. not intended to or capable of replacing human experts in the modeling process. No human-derived data was used.
Acknowledgments and disclosure of funding
This work is supported by Cancer Research and Roche. We want to thank Samuel Holt, Jonathan Crabbé, and anonymous reviewers for their useful comments and feedback on earlier versions of this work. We are grateful to Professor Yuanzhang Xiao for insightful discussions about the optimization algorithms.
10

---

References
[1] Ahmed M Alaa and Mihaela van der Schaar. Demystifying black-box models with symbolic metamodels. Advances in Neural Information Processing Systems, 32:1 11111134-1111 11314, 2019.
[2] MS Andersen, J Dahl, and Vandenherghe. CVXОРT: Python package for convex optimiza- tion, Version 1.1. 6, 2013, 2013.
[3] MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual. Version 9.0., 2019.
[4] L. Biggio, T. Bendinelli", A. Neitz, A. Lucchi, and G. Parascandolo. Neural Symbolic Regres- sion that Scales. In 38th International Conference Machine Learning, July 2021.
[5] Stephen A. Billings, Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Tempora Domuins. John Wiley & Sons, September 2013.
[6] Bongard and H. Lipson. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104(24):9943 9948, June 2007.
[7] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi. Sparse learning of stochastic dynamical equations. The Journal of Chemical Physics, 148(24))241723, June 2018.
[8] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
[9] Richard P Brent. Algorithms for minimization without derivatives. Courier Corporation, 2013.
[10] Steven L.. Brunton, Joshua I.. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932-3937 April 2016.
[11] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[12] 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. /927 [cx.NE], 2021.
[13] Jonathan Crabbé, Yao Zhang, William Zame, and Mihaela van der Schaar. Learning outside the Black-Box: The pursuit of interpretable models. arXiv:2011.0859 [cs], November 2020. arXiv: 201 1.08596.
[14] D.G Crighton. Model Equations of Nonlinear Acousties. Annual Review of Fluid Mechanics, 11(1):11-33, 1979.
[15] Carl De Boor. A practical guide to splines, volume 27. springer-verlag New York, 1978.
[16] Brian de Silva, Kathleen Champion, Markus Quade, Chan-Christoppe Loiseau, Kutz, and Steven Brunton. Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49):2 104, 2020.
[17] Bradley Efron, Trevor Hastie, lain Johnstone, and Robert Tibshirani. Least Angle Regression. The Annals of Statistics, 32(2):93, 2004.
[18] Lev D Elsgolc. Calculus of variations. Courier Corporation, 2012.
[19] F. G. (Friedrich Gerard) Friedlander. Introduction to the theory of distributions / F.G. Friedlan- der. Cambridge University Press, Cambridge, 1982.
[20] Samuel Holt, Zhaozhi Qian, and Mihaela van der Schaar. Deep generative symbolic regression.
In The Eleventh International Conference Learning Representations, 2022.
[21] Frank Hoppensteadt. An Age Dependent Epidemic Model. Journal of the Franklin Institute, 297(5):325-333 May 1974.

---

[22] James M. Hyman and Basil Nicolaenko. The Kuramoto Kuramoto-Sivaahinskky equation: A bridge between PDE'S and dynamical systems. Physica D: Nonlinear Phenomena, 18(1):113-126, January 1986.
[23] Kadierdan Kaheman, J. Nathan Kutz, and Steven L. Brunton. SINDy-PI: robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476(2244))20200279, October 2020. Publisher: Royal Society.
[24] Alan A. Kaptanoglu, Brian M. de Silva, Urban Fasel, Kadierdan Kaheman, Andy J. Gold- schmidt, Jared Callaham, Charles B. Delahunt, Zachary G. Nicolaou, Kathleen Champion, Jean-Chriss Saan-Chrittophe Loiseau, J. Nathan Kutz, and Steven L. Brunton. Pysindy: A comprehensive python package for robust sparse system identification. Journal of Open Source Software 7(69):3994, 2022.
[25] 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 Mathematicx, page 11, 2016.
[26] John R. Koza. Genetic programming: on the programming of computers by means of natural selection. Complex adaptive systems. MIT Press, Cambridge, Mass, 1992.
[27] Yoshiki Kuramoto. Instability and Turbulence of Wavefronts in Progress of Theoretical Physics, 63(6):1885-1905) June 1980.
Systems.
[28] 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, 3 31(11):40054442, 2020.
[29] Zichao Long, Yiping Lu, and Bin Dong. PDE-Net 2.0: Learning PDEs from data with a numeric- symbolic hybrid deep network. Journal of Computational Physics, 399:108925, December 2019.
[30] 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, pages 3208-3216 PMLR, July 2018 ISSN: 2640-3498.
[31] Lynn H. Loomis. Advanced calculus / Lynn H. Loomis and Shlomo Sternberg. Addison Weeleyy series in mathematics. Aidison-Wesledd Pub. Co., Reading, Mass. London, 1968. Publication Advanced calculus.
Title:
[32] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor. Model selection for dynamical systems via sparse regression and information criteria. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2204)22017000) August 2017.
[33] 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. Conference Name: IEEE Transactions on Molecular, Biological and Multi-Scale Communications.
[34] Mikhail Maslyaev, Alexander Hvatov, and Anna Kalyuzhnaya. Data-Driven Partial Derivative Equations 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, editors, Computational Science ICCS 2019, Lecture Notes in Computer Science, pages 635 641,, Cham, 2019. Springer International Publishing.
[35] Daniel A. Messenger and David M. Bortz. Weak SINDy for partial differential equations. Journal of Computational Physics, 443: 10525, October 2021.
[36] Daniel A. Messenger and David M. Bortz. Weak SINDy: Galerkin-Ban Based Data-Driven Model Selection. Multiscale Modeling & Simulation, 19(3):1474-1499 January 2021.
12

---

[37] Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondiej Čertík, Sergey B Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, and others. SymPy: symbolic computing in Python. PeerJ Computer Science, 3:e103, 2017. Publisher: PeerJ Inc.
[38] A. R. Mitchell and F. Griffiths. The finite difference method in partial differentiall equations. Wiley, Chichester [Eng.]: New York, 1980.
[39] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Counnapea,, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-leam: Machine Learning in Python. Journal of Machine Learning Research, 12:2825-288 2830. 2011.
[40] Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio Santiago. Soo Kim, and Joanne T Kim. Deep Symbolic Regression: Recovering Mathematical Expressions From Data via Risk-seeking Policy Gradients. ICLR 2021, 2021.
[41] J. Promrak, G. C. Wake, and C. Rattanakul. Predator-pre prey Model with Age Structure. The ANZIAM Journal, 59(2)1155-166, October 2017. Publisher: Cambridge University Press.
[42] Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-CODE: Discovering Closed- form ODEs from Observed Trajectories. The Tenth International Conference Learning Representations, 2022.
[43] Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125-141 March 2018.
[44] 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.
[45] Chengping Rao, Pu Ren, Yang Liu, and Hao Sun. Discovering Nonlinear PDEs from Scarce Data with Physics-Encoded Learning. ICLR 2022, 2022.
[46] Patrick A. K. Reinbold, Daniel R. Gurevich, and Roman O. Grigoriev. Using noisy or incomplete data to discover models of spatiotemporal dynamics, Physical Review E, 101(1)010003, January 2020.
[47] Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 364)e1600214, April 2017.
[48] 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):201604466 January 2017.
[49] Michael Schmidt and Hod Lipson. Distilling Free-Form Natural Laws from Experimental Data. Science, 324(5923):81 April 2009.
[50] Michael Spivak. Calculus On Manifolds: a Modern Approach To Classical Theorems Of Advanced Calculus. CRC Press, 2018. OCLC: 1029237047.
[51] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd. OSQP: operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4)6637-672, 2020.
[52] Trevor Stephens. gplearn: Genetic programming in python, with a scikit-learn inspired and compatible api, 2022.
[53] Robert Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Sta- tistical Society. Series B (Methodological), 58(1):267-288, 1996. Publisher: [Royal Statistical Society, Wiley].
[54] Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Al Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. 34th Confer- ence on Neural Information Processing Systems (NeurIPS 2020), 2020.
13

---

[55] Silviu-Marian Udrescu and Max Tegmark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 616))caay2661, April 2020.
[56] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C Carey, Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, lan Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, António H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261- 2020.
[57] Sachin Wani and Sarita Thakar. Crank-Nicolson Type Method for Burgers Equation. Interna- tional Journal of Applied Physics and Mathematics, 3:324-328, January 2013.
[58] G. F. Webb. Population Models Structured by Age, Size, and Spatial Position. J. M. Morel, Takens, B. Teissier, Pierre Magal, and Shigui Ruan, editors, Structured Population Models in Biology and Epidemiology, volume 1936, pages 1-49. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. Series Title: Lecture Notes Mathematics.
[59] Glenn F. Webb. Theory of Nonlinear Age-Dependent Population Dynamics. CRC Press, January 1985.
[60] Christopher K. Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
[61] Hao Xu. DL-PDE: Deep-Learning Based Data-Driven Discovery of Partial Differential Equa- tions from Discrete and Noisy Data. Communications in Computational Physics, 29(3):698-728, June 2021.
[62] Hao Xu, Haibin Chang, and Dongxiao Zhang. DLGA-PDE: Discovery of PDEs with incom- plete candidate library via combination of deep learning and genetic algorithm. Journal of Computational Physics, 418: 109584, October 2020.
14

---

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 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
[
set of numbers (1....,n)
N
set of natural numbers, i.e., {1, 2.........
No
a set of non-negative integers, i.c., {0, 1, 2, ... the dimension of the domain of vector field the dimension of the codomain of vector field denotes the smoothness of functions the maximum order of derivatives the size the dataset observed fields
M
N
K
D
S
the number of testing functions an open set in T the derivative of at multi-index, element of NM the order of a. a Σ M a;
Ω
i(t)
o:
a -order partial derivative with respect the 5th variable 0 0g "9.....0% ...0%
go
ck
set of functions with continuous partial derivatives 8 for all aa| K
E
an extended derivative, Definition |
an extended differential operator, Definition]
F
the functional used in the variational loss, Definition
(()))))
the variational loss, Equation sampling grid, Definition 10
u
true field, Definition an observed field, Defimiion dataset of observed trajectories the noise
D
Q
dictionary of non-degenerate extended derivatives vector describing the Ô-bound part of the VR-PDE
3
noise ratio
AL
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). closed-form expression is mathematical expression that consists of finite number of variables, constants, arithmetic operations, and certain well-known functions (e.g., logarithm, trigonometric functions). function called closed-form it can be represented by closed-form expression. E.g., f(z x² log(y) sin(3z).
15

---

Remark. In practice, we 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 4
Definition 5 (Muuli-indexx), An n--dimesiinaal multi- index a is an n-tuple
a = (a1,
an)
where Vi E [n] 0g No- Thus No- We define the order of a |a| = Σ Oz-
Definition 6. For any n--dimensional multi-index a we define a mixed derivative
8 0 "..".00
where 0.eh = 8 /0, is a th order partial derivative with respect 2; (the 7th independent variable). We call d non-trivial if [a] > 0. We denote the list of all non-trivial partial derivatives of up to order K as 0 k)
Definition 7 (Closed-form Partial Differential Equation). Let f be a closed-form real smooth function. We say that a vector field u -> R is governed by Wtth -order closed-form PDE described by f
f(z u(x), 00
(())) = Væ E Ω
(17)
where glk) u are all non-trivial mixed derivatives of uj (j [M] up to the h order. Definition 8 (Testing function). Support of a function Φ: $2 + R is defined as
supp { E Qt <<x) / 03
where B is the topological closure of B in Ω.
is called testing function it is CK function with compact support.
Definition 9 (True Field). We define true field on Ω as a vector valued function u : Ω + RV where each uy 12 R is C function.
Definition 10 (Observed field and sampling grid). We define sampling grid G to be finite subset of Q. Let u Ω 1 be a true field on Q. An observed field sampled from u on grid & is
function v G + RR of the form
(x) = uj (x) + <)( Va E GVj E [N]
where (2)) corresponds noise, realisation of zero-mean random variable.
Definition 11 (L1 sphere). Let n E N. We define n-dimensional L1 sphere to be subset of R=+1 defined as:
{ar E R2+1
||/||| 1} c R
(18)
Definition 12 (Standard simplex). Let N. We define standard nimmlee to be subset of R defined as:
#+1
{o E Ra+1
I = 11 2 201 [n 1}} C 1=+
(19)
B Variational-Ready PDEs
B.1 Variational Formulation of PDEs
In this section, we provide the standard variational formulation of PDEs for linear PDEs 11 Definition 13 (Linear differential operator). Let A be a finite set of multi-indices. linear differential operator L is defined as
L
Σ da
acA
where aa C is non-zero sufficiently smooth function of dependent variables. max |a then we call L an nhh order linear differential operator. If all o, are constants we say that L has constant coefficients.
The adjoint L,, denoted Lt, is a linear differential operator defined as
L Lu((x)
Σ
(-1)) gor (aa (x))()))
(20)
acA
16

---

Proposition 1 (Variational Formulation of PDEs for linear PDES), Let K E N. Consider scalar field u Q + R. such that CK2 Rthh order linear differential operator L, and a continuous function Ω - R. Let S2 - R be testing function. Then satisfies a linear PDE
L[u(x)] - g(z) -A EΩ
(21)
if and only if
J. [u(z)Lid())) g(x)))(=)) dx
(22)
for all testing functions Ø S2 + R.
Note that the integrals are always well-defined o 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 [8. We also need generalized version of the divergence theorem [31]
Lemma (Fundamental lemma of calculus of variations). Let K N. Q be open in R. and : Q + R be a continuous function. Then u is equal to 0 on the whole if and only if for all CK functions Ω - R with compact support.
Ja u(z)
Proof. If u is identically 0 on then all the integrals are trivially equal to O.
We now prove the converse.
Let as assume for contradiction that there exists a point 3y Q2 such that u(wo) A O. Without loss of generality we assume u u(zo) = c > 0. As is continuous there exists an open ball around zo of radius 8, denoted Be = fxr E Ω || so la < 6), such that Var E B u(x) > c/2 > 0.
Now let Φ be a CK function that is positive on BA, and elsewhere. Such a function can always be created by appropriately shifting and scaling ((x) e -110111=08) - [[lalaa<1)- Its support is
closed ball B6 so = (æe E Q
|| xo S 83 which is compact. Then
J
(()))(x))
u(=))(=)dz >
(23)
B2g
both u(x) and ((x) are positive B Thus found continuous function & with compact support such that:
[
u()))(())da #
(24)
Therefore u is identically 0 on Ω.
To make this section self-contained, we provide the statement of the generalized divergence theorem
31.
Theorem 2 (Divergence theorem). Let f2 be an open set in R& and f.g g.g be continuous on Ω = QU N2 and continuously differentiable Q. Then
S [[(((=))((=)d==== f( f(x) f(x)0 f p(f(z)g(z)dz
(25)
where v is normal unit vector to the boundary an.
In 1 1-dimensional setting, the statement of the theorem reduces to the integration by parts.
We can now prove Theorem [
Proof. Let us denote Ep- = (Otp: : upp hp). Then the PDE in Equation 9 can be written as:
P
Σ a)(2))00 (1 u(x))] g( u(x)) 0 Væ e Ω
(26)
17

---

The LHS is continuous as all ap and hp are smooth, is continuous, u E CK and |apl - K Vp E [P]. Thus we can use the fundamental lemma of calculus of variations say that the Equation 16 true if and only if
Jo
Σ
(2, u(x))] g u(x)) ((x)dx 0
(27)
for all testing functions Φ. We transform the LHS of Equation 27 using linearity to:
Σ ap(x)000 (h u()))])(()) J g( ((x))) p=1 n
(28)
Let us now focus on
/ (2(x0)000 (hpla, u(x))]]((x)
(29)
and let us denote ap (opl....,,0) (par) Then 00, 0001 0. and the expression can be written as
J 8 ... M [,,
(30)
Let us denote the support of d B. As is equal to zero outside of its support, can write the expression as
Je 8.50 1f.t....0 00 Dh
(31)
8
Without loss of generality, let assume that apl > 0. By the divergence theorem, this can be rewritten as
OM [h, (æ. u u(x)) (=))))][(,
(32)
B because the integral the boundary is equal to 0
J 11 ...0% [1,( (z)d))
(33)
as & has a compact support (and thus vanishes on the boundary). We can perform this operation 0pp times to shift the whole derivative 0Prlt to the second part of the equation and obtain
(-1)
0512 .......0a (hp u
(34)
Then we repeat this for other derivatives and end up with the following expression:
(-1)aa (-1)agm [ h,, u(x))00") ... (a (ap(z))(())dde
(35)
As the integrand is zero outside of B, this can be rewritten as:
(-1)lapl [ hp(2,, v(x))00 (a,(z)
(36)
or more compactly, using the functional defined in Definition as:
F(Ep: u. ②)
(37)
Therefore Equation 28 can be written as:
P F(E, u) () -000/11 (9(2,, u(z)))(()) dæ
(38)
p=1
Thus, we proved that Equation
26 is true and only if
P F(Ep, u, Φ) - / [0(2 (x M
(39)
for all testing functions
18

---

Table 5: Examples of equations which are Variational -Ready
Name
Equation
Linear
VR
Damped wave eq. with source
+ pu, KV - g(xz)
Gauss law
V = P/E0
Burger's equation
Wur Vulax = 0
X
Navier-Stokes equations
ut + (u. V)u - vV 1 - 1/VP + +9
X
Korteweg-De Vries equation
u Uxzz Giily = 0
Seramott--Svvshinsee equation
Ue + Uzz + Uxxxz + Wl2 = 0
Fisher's equation
U Kugx ru(l 2)
Liouville's 's equation
U + Uyy KEPu
Porous medium equation
U (²(1²) =0
Sine-Gordon equation
UST sin (u)
B.3 Examples
The examples of VR-PDEs can be found in Table
C D-CIPHER
C.1 Rewrite the inner optimization as a constrained least squares
Let us rewrite the objective in Equation
14
D
s
P
ΣΣ
F(BpEp
al4 / g(Ø, (z))
(40)
First, us observe that
F(SrEp (1)) .h) / 11ple,
= -pp2
(41)
we let z/ (d.0) E R defined as
(d,a) ( 1 al (1)) (0)
2p
(42)
5
Moreover, if we define (d.s)) RR
(dds)) / g(2, (z)))
(43)
we can rewrite expression 40 as
ΣΣ
Σ Spop ((d
(d(())))
(44)
d=1
Now, the sum over p can be written dot product between (d.s) E RR. and 3 E R We can also combine the sums over d and s. We obtain
Σ [
(g. 2(.) (d.e))
(45)
which is exactly the same as the objective in Equation 15
C.2 Pseudocode
The pscudocode of D-CIPHER is presented in Algorithm
19

---

Algorithm D-CIPHER
Input: Observed fields D { grid g
Input: Symbolic regression optimization algorithm O
Input: Smoothing algorithm S
Input: Testing functions {
Input: Dictionary Q p_ E (ap, Up:
Step
Output: Target PDE
a = =(v(d) vd E [D]
Step
initialize matrix Dx8 x R
2/d.3) - h, a(d) (x)) (x)) ga,
procedure Lgss(g)
initialize vector
R Dx5
w(d,s) + fg g(x à ( (z)dz
3 + COLLIE Z, w)
>> Section 7
L= 1 -
return L end procedure O(Los))
Step
initialize vector w Dx8
(d.s) t Jn (α, a (z)dx
B + COLLIE Z w)
>> Section
return 1p 3pEplu
u(x))
C.3 Testing functions
Testing functions should satisfy the following conditions:
1. Be sufficiently smooth (at least CK for K order PDE)
2. Compact support
3. Derivatives be computed analytically
4. Orthonormal
Conditions and follow directly from Definition 3 Condition 3 is necessary because we do not estimate the derivatives of the testing functions. Condition 4 follows from the result obtained by [42] that suggests that these functions should be a subset of an orthonormal basis of L2 space.
We use B-Splines 15 as the testing functions in our experiments because can control their smoothness and the derivatives are easy to compute. We scale and shift them appropriately so that they orthonormal.
Other testing functions are possible and examining them constitutes an interesting research direction. In particular, piecewise polynomials defined by (35] or various wavelets. Ideally, we would like to choose wavelets that form an orthonormal basis for the L.2 space, such as
Shannon wavelets smooth (Coo) but not compact
Meyer wavelets smooth (Co) but not compact (better rate of decay than Shannon)
Daubechies wavelets smooth (C for a specified K) and compact but they do not have closed-form expression.
Another interesting avenue of research would be to adapt the testing functions to the input data.
20

---

D CoLLie
D.1 Lagrangian
The problem that CoLLie is supposed to solve is a constrained least-squares optimization defined as:
minimize
||Az - bl12
(46)
subject to |||||
where E Rm has a full column rank, b E R., and x E IR" for some m, n E N. We consider the Lagrangian L : R" x IR + R associated with this problem [8] defined as T(2, = || - bl2 + >[1==11- 1)
(47)
Now let us define & + 12 as
2 arg min L(z, A) arg min ||A= b112 + \||=||1 zER
(48)
ZER"
The goal of our algorithm to find X E E such that ||
Let us define function
R + IR as
(()) - ||
(49)
The goal can be phrased as finding X E R such that g(>²) =
Let us note that 20) is just a solution to the ordinary least squares (OLS) problem with no constraints and its norm (()).
D.2 Extending LARS
Case 1. (()) 2 1.
If we assume that 2 is continuous then q is also continuous. From the continuity and the fact that (A) oo
and q())= 2 we infer that there exists a 2 such that g²)) 1. Moreover, for À 2 0) the problem in Equation 48 is the same as in LASSO 53
Therefore we just need to perform LASSO for different À and choose the one that gives the solution with L1 norm equal 1.
To do it in practice we use Least Angle Regression (LARS) 17 a popular algorithm used to minimize the LASSO objective. It generates complete solution paths, i.e., function e RR, - RR" defined as
c(A) = arg min A= bl + \||=||| ZER"
(50)
which is equivalent to for À 2 0. An illustration of LARS solution paths can be seen in Figure Each line corresponds to function which describes the coefficient for the ht covariate. The paths are defined from some An where all (Ap)) = 0 to À = 0 where c(0) = 2()) In other words, the solution paths cover the whole range of constraints from the strictest, effectively imposing the L1 norm of z to be O, up to 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 (Ao > ... > An = 0) where the slopes change. We calculate the norm at each of these points, ||c(Ax]lı, and find j E [n] such that ||c())-1)llı a linear function on [ Ay-1] and we know both c(Ay-1) and c(2y), we can effectively search for A Ay-] such that |c(A)] = 1. The search can be performed by any root-finding algorithm.
S [[c(A))]l1- As each c,
We use Brent's method 9
Case 2. 0 < (() 1
This is much more difficult as it corresponds solving the problem in Equation 48 for À < 0. The solutions given by the LARS algorithm are too small. In fact, the solution with the biggest norm is c(0)= (0), the OLS solution, with norm exactly ((0)= < 1.
To address this challenge, we propose the following heuristic. We extend the solution paths generated by LARS beyond X 0 for À < O. We assume that the paths will continue to be piccewise linear and that they will keep the slope they have in the last interval [An = O, An-1]. Let us denote this slope as
Aç (0) Cs (Aa-1)
(51)
0 An-1
21

---

This is graphically represented in Figure 6 Formally, these extended paths, E: IR - R are defined as:
[c((A), (A)) = (c(0) + AAc;
A X 0
(52)
À < 0
Now, we want to find À 0 such that |[(A)]] To achieve this in practice, 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 ( is on the wrong side e, (A)Ac, > In other words, being on the wrong side just means that the path yet to cross the x-axis if we keep decreasing A. We can easily find X such that for all X < X all a (A) are on the right side (none of the paths will ever cross the x-axis).
(0 c,(0) {
lie [n] A (0)40 > o0}
X min
(53)
Ac
If ||(\) 2 we just need to search the interval [X',0) for such that ||ε(\\]] 1.
If (e())la < 1 then we need to search X < X. However, by definition, for all < X, all (A)
the right side. That means ||ε()]] as function of is just linear function the interval (-o0, X'). To see that, let us observe that
|(())) = Σ [c(A)] = Σ (i))(((
(54)
Additionally, for < X' all (A() are on the right side, so we have sign((a(A)) = sign(Ars). We rewrite ||(A)]]ı as:
|[(()) = Σ (G)gg((A))((())) (Aci))
i=1
sigm(Ar:)cs))
Σ sign( Dc:)) Ari
(55)
sigm(Ac))
Σ
Aczl x
Therefore the solution can be found using the following equation
x = X' +
(56)
Σm=
Acj
Case 3. q(0) = 0. In that case, we just return a precomputed solution to the problem
minimize
||
(57)
subject to | 11
1
which we compute by subdividing the problem into 2n quadratic programs and solving each of them separately using CVXOPT algorithm (2] as described in Section 7
D.3
Comparison
We perform a comparison between CoLLie and algorithm based on CVXOPT as described in Section 7 CoLLie uses different procedures depending on the problem, thus for fair comparison we generate an equal number of tests falling under Case and Case (two main cases) as described in [.2] To achieve that, we generate random m X n matrix A where each entry is sampled from
standard normal distribution. We use m = 1000 and n ranging from to We then generate a vector 1 E R such that each entry is sampled from uniform distribution on [-0.5, 0.5). Then 2 is normalized have II norm equal to 1. We sample number / from a uniform distribution on [- 1,3)
and multiply & by to obtain z' 12. By this procedure, are guaranteed that cases || || <
22

---

(A)
(B)
Coefficients centtts
Coefficients
An
A1
An
An
An
Ax
An
Aq
An
Figure 6: Panel A shows and example solution paths calculated by the LARS algorithm. Panel shows their extended versions defined in Case in D.2 The x-axis is reversed, decreases moves to the right.
and ||=1|| 2 are equally likely. We then let Az'. The task then to find with norm equal to that minimizes |Az bl.. As g(0) = 1|>1|11 I, Case I and Case 2 are equally likely.
We perform 1000 experiments for each n. As losses for optimal solutions can be widely different scales, we report the relative етrоr 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 on a single computer with an Intel Core i5-6500 CPU (4 cores) and 16GB of RAM.
E
Experiments
E.1
Ablated D-CIPHER
The ablated version uses the standard MSE loss with estimated derivatives and thus solves the following optimization problem:
min min
ΣΣ Σ @pEplol())(()) g(x, v100 (x)))
(58)
|| d=1 FEG
where (d) is the observed field, G is the sampling grid, and Eg (xx)()((x) requires derivative estimation.
The ablated version uses the same symbolic regression algorithm to search over closed-form y and CoLLie for the inner optimization.
E.2 Implementation
Step 1. For the homogeneous heat equation and Burgers' equation we use dictionary Q {u, d.u, O,u, O2 (u²), Su, (u²)). For Kuramoto-Sivashinsk equation {u, Ōzu, O,u, , 0, (u²), oe, u, a (u²), Hú, 0z(u²) 2) Ozu, & (u²)). For the damped and forced har- use the dictionary {0,, 0f). and for the wave and heat equations, we use
dictionary Q
monic oscillator {Ōc. 0z, 02, 0,02, 02).
Step 2. Field estimation is performed using the Gaussian Process Regression from the Python library scikit-learn 39. The kernel chosen to the RBF kernel (60) with an added White kernel to account for noise. The observed field is initially standardized by subtracting the mean and dividing by the standard deviation. Then the CaussianProesRRRgessoo is fitted the data. The estimated fields are generated by predicting the values of trained Gaussian Process full integration grid
23

---

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 gpleam [52]. We custom fitness function that solves the inner optimization problem in Equation| [4] This inner optimization is performed by CoLLie (Section The integration is performed using Riemann sums.
Ablated version 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 estimate the derivative in one of the coordinates for all points in the sampling grid. To obtain higher-order derivatives, a Gaussian process is fitted again and the derivative is again calculated using the finite difference (possibly in 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.., (le - 5, 1e5).
Gllearn. We do not perform parameter tuning for the gplearn library and use the same parameters as in D-CODE 42 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 [52] 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 different experiments we record the loss on widely different scales. To prevent tuning this parameter for every experimental setting introduce penalty that work different scales. The parsimony coefficient is chosen manually by performing experiments for a few values. The value used in the experiments is 0.05.
The set of allowed mathematical operations is: {+, x, sin, exp, log)
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 10 testing functions and the integration step 0.01. For the inhomogeneous heat and wave equations we use 100 testing functions and integrate on a 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 as described above. The interval used in the finite difference method estimate the derivative was chosen to be: 10-3
24

---

E.4
Choice of equations
Equations used in Section 8.1 are canonical equations from physics that often appear in other works about PDE discovery [47] The homogeneous heat equation is a second-order PDE that models how heat diffuses through region. It contains the dissipative term 02u. Burgers' equation is second-order PDE used, for instance, in fluid mechanics or nonlinear acoustics [14]. It contains the advection term udgu and the diffusion term that prevents shock formation. Surramoto-Sivashine a fourth-order PDE used in modelling reaction-diffusion systems 27 and is known for its chaotic behavior
equation is
22
In Section 8.2 we chose equations of physical significance that have an interesting Ò-free part that is not a linear combination (as discussed in Section 3). That makes them impossible to discover by the current methods. A forced and damped harmonic oscillator is a second-orded ODE. Although D- CODE 42] can discover any closed-form first-order ODE, it cannot be used to discover second-order ODEs. Thus there is currently no algorithm capable of discovering this equation. Inhomogeneous heat and wave equations are second-orde 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 a term. Thus even without the source term, most of the current methods cannot be applied directly to discover this equation.
E,5 Data generation
Homogeneous heat equation. The fields were generated by solving the equation Ō,u - 0,00,u 0
(0, = 0.25) with Neumann boundary conditions θ, ut, 0) = (t, X) 0 and an initial condition x0) = uo (z), where Un is randomly sampled from a Gaussian process. The equation is solved using the implicit BTCS scheme [23] with steps öt 0.001 and óx = 0.001. The observed field is generated by sampling (t, x) E [0. T] x X). evaluating the true field u(t, x) and adding Gaussian noise. T and X are used in the experiments.
Burger's equation. The fields are computed by solving d.u + u0,u - 0, 8 = 0 (0) = 0.2) with an initial condition u(0,xz) uo (z), where uh is randomly sampled from a Gaussian process, and with Dirichlet boundary conditions u(t, 0) = un (0). utt, X) = uo (X). The equation is solved the using Crank-Nicolsor scheme (57] with steps òt 0.002 and óx = 0.002. The observed field is generated by sampling (t, x)) E [0, 7 x [0, X]. evaluating the field u(t,, x) and adding Gaussian noise.
T = 2 and X = 2 are used in the experiments.
Kuramoto-Sivashinsky equation. The solution the same as the one used in (47. The observed field is generated by sampling (4, x)) [0, T" x [0, X). evaluating the field w(l, 2)) and adding Gaussian noise. 7 100 and X = 100 are used in the experiments.
Damped and forced harmonic oscillator. The true fields are created by analytically solving the equation of (t) 20200 (z), + 0lg(()) = 0g sin(), where 01 0.5, 02 4.0, 03 5.0, 0g 3.0, with random initial conditions for u0) and (,)( The observed fields are then created by sampling E [0, 27. evaluating u(t), and adding Gaussian noise. = 2 was used in the experiments.
Inhomogeneous heat equation. The true fields are computed by solving oeu(t, x) - (t, () Oze where 0, = 0.25, 0₂ = 1.25, 0g 1.8, with Neumann boundary conditions (t, 0) O.,((, X)= = 0 and an initial condition w(0, x) = ux))). where to is randomly sampled from Gaussian process. The equation is solved using the implicit BTCS scheme (25) with steps ôt = 0.001
and dz 0.001. The observed field is generated by sampling (t, x) [0, T] [0, X), evaluating the
true field xtt, and adding Gaussian noise. T = 2 and X = 2 are used in the experiments.
Inhomogeneous wave equation. The true fields are computed by solving (fu(t,z) - 0 (, Bzet (fst), where 81 1.0, 0₂ = 2.0, 0g 3.0, with Dirichlet boundary conditions utt, (0), u(t, x) = Mx)X), where ue is randomly sampled from Gaussian process and specifies the initial condition u(0, ) = =n(x). The equation is solved using the Implicit Difference Method (38 with steps 8t 0.001 and 5ι = 0.001 The observed field was generated by sampling
(I, x) E T) x [0, X). evaluating the true field u(t, x:) and adding Gaussian noise. T 2 and
X were used in the experiments.
25

---

E.6 Experimental settings
Homogeneous heat equation
Noise ratio (Øg): 0.001, 0.01, 0.1
Number of samples (D): 10
Domain (It): [0,2) [0.2]
Grid (9): {0, 0.07......... 2} x (0,0.07,
Burgers' equation
Noise ratio (og): 0.001, 0.01, 0.1
Number of samples (D): 10
Domain (S2): [0, 2] x [0.2]
Grid (9):
2} (0,0.1 (0.......... 2]
Kuramoto- Sivassik-S equation
Noise ratio (ox): 0.001, 0.01, 0.1
Number of samples (D): 1
Domain (St): [0, 100) [0.100]
Grid (9): {0, 1.67.... 100} x (0, 1.67, 100}
Damped and forced harmonic oscillator Default values are in bold.
Noise ratio Ro 0.001, 0.005, 0.01, 0.1, 0.2, 0.5 Number of samples (D): 1, 2, 5, 10, 15
Domain (R): [0,2]
Grid (9): {0, 0.08, .........)... Inhomogeneous heat equation
................. 2).
(0. 2). (0. (0,0.2.. ... 2].
(0,004
Noise ratio (σ 0.05, 0.1, 0.2
Number of samples (D): 10
Domain (I)) [0.2] 2] [0.2]
Grid (9): {0, 0.07. 2} x 00,0.7......
Inhomogeneous wave equation
Noise ratio (og): 0.001, 0.01, 0.015
Number of samples (D): 10
Domain (?)): [0, 2] [0.2]
Grid (G): f0, 0.07. ..... 2} x ..,
E.7 Benchmarks
We use the implementation PDE-FIND and WSINDY from the PySINDY library (24 [6]
An important hyperparameter in PDE-FIND and WSINDy is the library + used. We impose that Q and e U {ō,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
θ {u, Oz, 8. u, ud ud.u).
26

---

For the Kuramoto- Sarramott--ivsssissky equation, we use a library
@ (u, Oz u. Bdu, ozu, ud,u, udLu, uO u0z).
In the experiments, we have not optimized for the derivative-free part as identically equal 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 [42] considered their functional forms, i.c. expressions where all numeric constants are replaced by placeholders. By this measure, functions sin(3x) and sin(3.52) match as they have the same functional form (Cr), where C' is a placeholder.
However, this definition is quite restrictive because functions sin(3.r) sin(32) 0.0001, 1.001 sin(3z), and sin 0.001) all have different functional forms.
We consider an open challenge to design a good metric that would meaningfully reflect whether the correct equation is discovered. We propose the following.
For target function f. we consider its augmented form f. defined f(x) Cif(Cza + C4) + C2. where C, are placeholders. Then all numeric constants are turned into placeholders as well. the end, we combine the constants. For instance, C1 + C2 becomes just C3.
As an example, let us consider a function f(z) 1. The augmented functional form is created in the following way:
1. Augment: C1 L.3e + C2
2. Replace: C x C3eCo((Cxxxtt)) + C2
3. Combine: Cie Cg2+Ca C2
We perform this procedure for the target function. We can now take the standard functional form the candidate function and check whether it matches the augmented functional form of the target function, taking into account that some the constants might not be present in the candidate expression.
To aid in this procedure, we use a Python library for symbolic mathematics, SymPy [37
E.9
Computation time
The average computation time for a single experiment with the damped and forced harmonic oscillator 281 seconds with standard error of 4.5 seconds. The average computation time for a single experiment with an inhomogeneous heat equation is 68 minutes with standard error of 38 seconds. This time measured on single computer with an Intel Core i5-6500 CPU (4 cores) and 16GB of RAM.
The experiments
simultaneously computers like the described above. The total time for all experiments (all seeds, all equations, all experimental settings, and both versions of D-CIPHER) is 65 hours.
E.10 Licenses
The licenses of the software used in this work are presented in Table 7
Additional experiments and discussion
F.1
Sharp--lolkaaMeeeediiik model
We test D-CIPHER on population model called the Shamp-Looke-Meeendriikk model [59 The model is described by the following equation
azu(t, a) d.u(t, a) m(a)u(t,a) =
(59)
27

---

Table 7: Software used and their licenses
Software
License
gpleam
BSD 3-Clause "New" "Revised" License
evxopt
GNU General Public License
evxpy
Apache License
sympy
New BSD License
scikit-learn
BSD 3-Clause "New" "Revised" License
numpy
liberal BSD license
pandas
BSD 3-Clause "New" "Revised" License
scipy
liberal BSD license
python
Zero-Clause BSD license
pysindy
MIT License
where m(a) is age-specific mortality rate. We choose m(a) = 2e We set = 1.5 in the experiments.
We note that the derivative-free part of the target PDE cannot be expressed as a linear combination 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 mortality rate that can be represented as closed-form expression. We show comparison between D-CIPHER and the Ablated D-CIPHER in Table
Table 8: Simulation results for the Laare-Lotta-MMKenddee model. We report the success probabil- ity of discovering the 0-free part and the Average RMSE of the 0-bound part. Standard deviations are shown brackets
Sucess Probability
Average RMSE)
Method
0R 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)
F22 Discovering systems
Discovering systems PDEs is much harder problem than discovering a single PDE. One of the issues is the fact that would call indeterminism. It follows from the following fact. If vector field u is a solution to two differential equations fi and fz. i.c.,
(a u (d) (æ), a al u (d) (m)) 0 Va: E Ω fa(æ, u (x), kl (1) (x))
(60)
Væ E Ω
then it is also solution to any linear combination of these equations, i...
A1 fi (x,, u ( (az), (n()
(z:) X2 x f2(æ, u (d) (zz), 0 k) (d (x)) Va:
(61)
for any Å1 A2 R. Moreover, equations can sometimes be differentiated to yield equations. Let us take as an example the Cauchy- Riemman equations defined as:
0241 - 0y = 00 + 0y =
(62)
Let assume that we have true vector field (μ, Uz) that satisfies both equations. following equations are also satisfied
Then the
0_40 0,u2 0,42 + 0,u1 =
(63)
0,01 - 0,u2 0242 - 0,u1 = 0
We can also differentiate the Cauchy-Riemann equations to arrive at the Laplacian equations for u1 and u2.
8.41 0fu1 0
(64)
822 02 U12 0
28

---

We could also combine the first-order equations with second-ord order equations or consider even higher- order derivatives. Although all these equations are compatible with our vector field (u,, 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 encoded to guide the algorithm into the right equations.
Current methods do consider systems of equations or consider system of equations of very particular form. In the latter case, each equation models a derivative with respect to time of a different scalar field. The system is assumed to look like this:
u, = fi u(x)) klk u(x))
0,u2 = fa u(x), a u(=))
(65)
0uu fz u klki u(x))
In addition, the LHS of these equations is often assumed to only contain spatial derivatives.
D-CIPHER can be used to discover some systems of differential equations if enough prior knowledge 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.
Finstly, we note that D-CODE 2 has been shown to discover system of equations that looks like Equations 63 when all equations are first-order ODEs. As D-CIPHER reduces to D-CODE when applied to first-order ODEs, it is also capable of discovering such a system.
We demonstrate that D-CIPHER able to discover both Cauchy- Riemann equations if we use two different dictionaries. Each of the dictionaries yields differentt equation. Additionally, it is well-known fact that if a vector field satisfies Cauchy- -iich-n equations then the constituent scalar fields harmonic, i.e., they satisfy Laplace's equation. Based the same dataset are also able to discover both Laplace's s equations given another set of two different dictionaries. We note that Laplace's 's equation does not contain Ô₂ term, 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 0y = 0 we use Q {82u1, u, Pun 0,
For 0,u2 0y O use Q2 {θ, By u1 u2 07
For 0.u1 0u 0 use Q3 {0. 0 u1 02 0f }
For Lug- Hu 0 we use Q 10, u2 02u2, 86
In the experiments, we have not optimized for the derivative-free part it identically equal 0 in all equations.
F.3 Increasing the size of the dictionary 2
We test D-CIPHER on the Sharp--Loka-MeKenndiie
model (Equation 59p with different dictionaries. We start with a small dictionary O1 {Ō,u, 0,u), ), and we create every new dictionary from the previous one by adding one more extended derivative. The final dictionary contains 10 elements, (u²). (,(22) OIu = {0,u, u, 0,u, S.u, 0.u, 0,0,u, 0g (u²), 0
0 (u²), (u²), 0z(2²) (u²))) The Average RMSE of the O-bound part is shown in Figure 8 We do not observe any increase in average error. Note that in these experiments, we just focus on the 0-bound part and do not optimize the Ô-free part. The relationship between the computation time and the size of the dictionary represented in Figure
9
F4 Computational complexity
We want to emphasize that PDE discovery not a 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

---

0,41 8yH2
2x4e + 8,41
0.06
0.06
Avg. RMSE
0.04
Aug. RMSE
0.04
0.02
0.02
10
10-2
16-1
0.0
10-3
10
16-1
Noise
(ow)
aºuı 8. 0
8² + 8,
0.45
0.24
Avg. RMSE
Avg. RMSE
0.16
0.3
0.15
0.0
10
10-2
10-1
0.0
10-3
10-7
10-1
Noise ratio (ow)
Noise ratio (ox)
Figure 7: Simulation results for both Cauchy-Riemant Riemann equations and two Laplace's equations. We report the Average RMSE of the 0-bound part in different noise settings
0.015
Avg. RMSE
0.010
D-CIPHER
0.005
0
2
4
10
Size of the dictionary Q
Figure 8: Simulation results for Sharpe-Lotka model. We report the Average RMSE of the 0-bound part for different sizes of the dictionary Q
Computation in D-CIPHER is performed Step 2 and Step
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²3) where n is the number of data points in a grid G. D-CIPHER is specifically have not encountered major computational issues
designed to work for sparse and noisy data, while performing Gaussian process regression. We also note that significant progress has been made in adapting Gaussian processes for datasets with many data points [28].
Step 3. D-CIPHER consists of two optimization loops. The outer optimization performed by symbolic regression algorithm (in our case genetic programming). The inner optimization is performed by CoLLie (Section 7). Searching through a space of closed-form expression requires testing many candidate equations. D-CIPHER is designed to work with many different algorithms for symbolic regression and it is advised choose an algorithm that can search through this space most efficiently.
30

---

0.28
2.012
Compution Inputtttonnno time
0.30
0.08
0.24
2
of the dictonnrry
Figure 9: Computation time of the experiments in Appendix E.3 dictionary Q influences the computation time
The plot shows how the size of the
CoLLie was specifically designed to solve the optimization problem quickly possible with minor accuracy trade-off (see Figure 2 It is based on LARS which has time complexity O(mn²) 17 where the number of samples and n. is the number of features. In case, n 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 10. Overall LARS is performed very quickly. The additional steps in CoLLie require only a few arithmetical operations (linear in n) and a possible root searching of single variable function that efficiently implemented using Brentq algorithm
Other important parts of the algorithm are the numerical integrations. One such integration is performed at the beginning of the algorithm to compute the matrix Z (see Algorithm []. It does not contribute much to the computation time as is performed only once. The other integration is performed for each candidate equation compute vector Also, substantial time is spent computing the values of the candidate function y used the integration. Fortunately, both of these operations be implemented as vectorized operations which are designed modern hardware. We also discourage overly long equations for g (check the discussion in Appendix E.3p to limit the number of operations performed.
very efficiently
F55
Challenges of derivative estimation
One of the advantages of D-CIPHER compared to other methods is the use of the variational formulation of PDEs that allows it circumvent derivative estimation. This is important as derivative estimation is challenging, especially in noisy settings with infrequent sampling. The problem becomes more pronounced the higher the order of the derivative. To demonstrate these issues, we perform series of synthetic experiments.
Qualitative study. First 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 10 (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 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.
that the estimated derivatives do not match the
Quantitative study. We investigate this relation quantitatively for the damped and forced harmonic oscillator and the wave equation. For the oscillator we generate 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
31

---

0.5 (8)
2.
Observed
Ground trith Finlt D fffrrencc
0.0
0.2
-1.0
5p ine
1.5
GP
0.4
0.0
1.0
0.2
1.0
0.0 1.0
X
X
Figure 10: Panel A shows that estimation of the true trajectory can 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.
mean squared error. The results are shown in Figure 11 (Panel A). We can see that the error increases the higher the order of the derivative. For the equation, perform similar experiment but this time we estimate different mixed derivatives. We consider any mixed derivative 0f where i (0, 1. We demonstrate the results in Figure 11 (Panel B). We observe that the error increases
the higher the order of the derivative (i ++++) j).
(A)
(B) a
ax
a
Finite Difference
250
GP
0.01
200
RMSE
150
t
0.05 0.5=0..00.220...0 1.643..73
100
50
a
0
1
2
3
Derivative order
Figure 11: Panel demonstrates the error between the estimated derivative and the ground truth increases with the order of the derivative (performed for the damped and forced harmonic oscillator). Panel B shows that the same happens for the equation. The (i.j) entry of the heatmap should be interpreted as the RMSE between the estimated derivative 0 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 we summarize some of them and describe how our work addresses them.
F.7 Significance of the new notions
In this section, want to justify that the new notions we introduce (evolution assumption, linear combination form, derivative-boumt part, derivative-fee part, Variational Ready PDEs) are important theoretical contributions that help understand the landscape of different PDEs from the machine learning perspective.
-eaadd
The definitions we introduce let us characterize different classes of PDEs. These new notions complement the standard recognized PDE classes such as semi-linear, quasilinear, hyperbolic, etc. These standard classes were introduced predominantly to characterize the solving techniques or the properties of the solutions, whereas the notions we introduce relate to the difficulty of discovering such equations from data.
32

---

Table 9: Some challenges of PDE discovery and how we address them
Challenge
Novel contribution (how we address the challenge)
Noisy measurements
We use the variational loss function (Equation 11
Large diversity of PDEs
We allow any closed-form Ö-free part and require no evolu- tion assumption (Section 3
Leariing compact equations
We penalize long O-frre parts and use LI normalization for the 0-bound part (Section 6 Appendix H.3
Efficient search
We develop a quick algorithm, CoILiee for the inner opti-
mization (Section 7
Chaotic systems
We do not estimate the initial conditions [ or perform for- ward time stepping [29] which are computationally unstable for chaotic systems
Significance of linear combination form and the evolution assumption. In Table 10 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 to 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 13
Significance of O-bound part, O-free part and Variational tive estimation has been one of the main challenges of PDE discovery. The variational formulation allows to circumvent derivative estimation and thus is more robust to 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 part and the derivative
PDEs. The difficulty of deriva-
the derivative- variational formulation. That allows to define Variational-Reade PDEs as currently the broadest class of PDEs that admit variational formulation. We believe methods requiring derivative estimation underperform in settings with high noise. This definition outlines current limits of any method that circumvents derivative estimation that way.
part of the PDE due to their significance for the
an important contribution as
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.
LC EA Equation form
Optimization problem
Examples
Yes Yes 0.u, Ep Opfo(z, u(x), a (((())
Relatively casy. Can be for- mulated as finding sparse so- lution to linear least squares (similar ridge regression) Medium difficulty. Can be formulated P separate linear least squares prob- lems as above Medium difficulty. Find using symbolic regression Very difficult. Requires other assumptions (in this work, Equation 13)
10
Yes No En Opfole, u(x), alk) =0
F3
No Yes Oeu, g(x, uIx))
No No f(x, u(x), a ( u(x))
D-CIPHER
F88 Comparison with D-CODE
To clarify our contributions, we compare D-CIPHER to D-CODE 42 which we think closest in spirit to the algorithm we developed as both of them allow any closed-form derivative-free part and
33

---

use variational formulation to circumvent derivative estimation. We also note how the new notions we introduce help us compare the two works.
Algorithm presented in D-CODE 42 can discover any first-order explicit closed-form ODE, i.e., an equation of the form.
(zur(t) = g(
where 9 closed-form function. uses the variational formulation of ODEs to circumvent derivative estimation.
Table 11: Comparison between D-CIPHER and D-CODE 42
Property
D-CODE 42
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 Any closed-form function No
Learned: Epe 3pE
Derivative-free part
Any closed-form function No
Derivative estimation
D-CODE 42 can be considered special case of D-CIPHER. We can recover from D-CIPHER by choosing the dictionary to contain only one element, i.e., Q = 0,ug.
As the derivative-bound part is fixed, every ODE of that form admits the variational formulation. This not true for PDEs as there are derivative-bound parts that might prohibit the variational formulation. Thus D-CIPHER required careful consideration of the appropriate class of equations search over.
As D-CIPHER needs to find both the Ò-free part (function g) and the 0-bound part, the optimization problem is much more complicated. That is why we restrict the derivative-bound part of the PDE to be spanned by terms from the pre-specified dictionary and develop an efficient optimization algorithm, CoLLie. We emphasize that, is the for [42], do not put any constraints on the derivative-fres part of the PDE, apart from it 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 significantly more difficult than the one considered other works. The space of PDEs we consider is much more complex than the space of PDEs in a linear combination form. In other works (e.g., [47), (35,, the PDE is basically a vector P and the discovery task mostly reduced to finding a 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.P, CFF(M N), where CFF N) is a space of closed-form expressions in M N variables. This space combinatorial the functional form and continuous in real constants. This makes it very challenging derive any error bounds.
Homogeneous heat equation
Burgers' equation
Kuramoto-Sivane -irassio-Sivashinsky equation
0.5
0.25
0.25
80 RASE
0.4
Peg. RAGE
0.2
0.15
0.3
0.15
Nearest
0.2
0.1
Linear
0.05
0.1
0.05
Cubic
10-3
10--
0.0
10
0.0
10-2
10
(du)
Noise (dt)
(as)
Figure 12: 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.
34

---

E.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 [60]. this section, investigate other common interpolation algorithms. We implement D-CIPHER with different estimation algorithms. In particular, we compare Gaussian Process (GP) against: Nearest point interpolation (Nearest), Linear interpolation, and Cubie interpolation (Cubic). The implementation details of these three algorithms can be found in the scipy (56) documentation. The results are presented in Figure 12 We that the estimation algorithms that produce smoother functions (GP, Cubic) tend to give better results.
E.11 How comprehensive does the dictionary Q need be?
Table 12 shows 17 different differential equations and what terms they need in dic- tionary. All of them can be discovered with a dictionary of size 10, Q = {8,u, u, 0,0,u, 0, 0,02, 02u, U,u (u²²,, 0, (u²)), and 9 of those equations can be de- scribed with a dictionary containing just 4 terms, Q = {dzu, Ōzu, 8fu, P.u). It is thus likely that such small dictionaries equations have similar derivative-bout bound parts and differ by derivative-free parts. Being able to discover any closed-form derivative-free part is what makes D-CIPHER stand out.
sufficient discover most of the well-known equations. Many differential
Table 12:
This table lists various Variational- Reaatd PDEs (some of which are discussed in the paper) and shows which extended derivative terms they require. All of them be discovered with dictionary of size of 10, and 9 of them can be described with a dictionary containing just 4 terms. thus likely that such small dictionaries are sufficient discover most of the well-known equations.
Equation
8.u. u
U.u
0,8,
0f u
ULu u
8.u
8g (w²)
(u
Heat equation
Wave equation
Burger's equation
SLM model
Damped harmonic oscillator
Suramoto-Sivashing equation
Benjamin- Babjamin-Bonn-Maaee Bona-Mahony equation
Boussinesy equation
Chafce-Infante equation
Damped wave equation
Fisher's equation
Hunter-Saxton equation
Klein- Gordon equation
Korteweg- -Deeeg--eeeee Vries equation
Liouville's equation
Sine-Gor Sine-Gordon equation
Sinh-Goree equation
F.12 How to take advantage of the observed derivatives?
D-CIPHER can make use of the observed derivatives (if they are available) by adapting the dictionary. Consider a setting with dictionary Q {Ō,u, 0_u, a,O_u, uuu, oe u). If happen to have the measurements Ú,u then we can introduce a new variable e = iuu and change the dictionary to Q = (u. u.u, 0,v, d.u, Lu). Note, have performed experiments where the dictionary contains more than one dependent variable in Appendix|F] F.2 With observed derivatives, we can also enlarge the space of Variational-Read PDEs by allowing g (0-free part) depend
well.
F.13 Impact on real-world problems
D-CIPHER is especially useful in discovering governing equations for systems with more than one independent variable. For instance, spatiotemporal data temporal data structured by age or size.
35

---

In particular, we envision D-CIPHER to be useful in modeling spatiotemporal physical systems, population models, and epidemiological models.
Spatiotemporal physical systems. D-CIPHER may prove useful in discovering equations governing the oceans or the atmosphere. For instance, some places actively add remove CO2 from the atmosphere. These "sources" and "sinks" are likely to be described by a Ô-free part which D-CIPHER specially equipped discover. Similarly, with the ocean temperature where O-free part can describe a heat source. Another area of application can be modeling seismic waves across the earth ' crust. Here the O-free part can describe the vibration source (e.g., an earthquake).
Population models. Population models can be used in agriculture to determine the harvest or for pest control to predict their impact the crop. They have also been used in environmental conservation to model the population of endangered species. Population models have also been used in modeling the growth of cells to better understand tumor growth. Moreover, understanding the evolution of population pyramid for specific country may prove invaluable in ensuring its economic stability. As in all these scenarios, the rates of growth and mortality are likely to described by Ô-free part, D-CIPHER is uniquely positioned to discover such equations as an aid to human experts.
Epidemiological models. Epidemiological models are crucial during pandemic for better planning and interventions. For many diseases, the rates of mortality and infection are age-dependent. Thus, modeling the spread of disease using PDEs (rather than ODEs) might provide superior results.
E.14 D-CIPHER in practice
Below we discuss the things to consider while using D-CIPHER.
The order of the differential equation. One of the first considerations should be to choose the order of the differential equation K. For many dynamical systems, K 2 sufficient unless we expect very complicated behavior. Then, considering K = 3 or even K = 4 may be warranted. Note, that we show that D-CIPHER can discover a fourth-order PDE (Kuramoto- -iinontt---vssiiisyy equation) in Section
Homogeneous equations. Before searching through the whole space of closed-form, g (derivative- free parts), we can consider whether the equation we want to discover may be homogeneous. These experiments the restricted search space can provide quick insights before searching through all closed-form derivative-free parts.
Terms in dictionary.
For given order of a differential equation K, it is a good idea to include all standard differential operators up to order acting on all the variables. For instance, for K = and M 1 could choose PDEs with constant coefficients up to that order. To allow for non-linear PDEs we can include a term like (u²) 2u0,u u that often describes advection (as in Burger's 's equation).
{a,u, u,u, a,0, u, 0zu, 0fu). That allows to cover all linear
Dictionary steering when dealing with many dependent variables. When we deal with a system of PDEs rather than single PDE choosing a dictionary is increasingly important. As we explain in Appendix F.2] discovering whole systems of PDEs is very challenging and D-CIPHER is not designed to do out of the box. However, we show how that can be done in certain situations. We can steer what kind of equations are discovered by choosing the terms the dictionary.
Estimation algorithm. Estimation algorithms make different assumptions on the data- generating process and should be chosen based on domain expertise. As we show in Appendix F.10] algorithms that produce smoother functions, such as Gaussian Process regression and cubic spline interpolation, tend have good results. We can consider the advantages and disadvantages of these methods. For instance, Gaussian Process regression works very well for smooth signals. However, it is computationally intensive and might not perform well the signal not smooth enough (it has abrupt changes). Spline interpolation, on the other hand, is faster and more appropriate for less smooth signals, but it might introduce certain unwanted artifacts because of using cubic polynomials to interpolate the data.
36


Reviews:
[{"review_id": "FIxl01zRV2", "reviewer": "Reviewer_kC9e", "summary": "The paper proposes a new way of discovering closed-form Partial Differential Equations (PDEs) from data. This especially aims at high-order PDEs, especially when the specific form is not pre-assumed and there is a lack of observations on derivatives. The key idea is to represent the unknown PDE with terms that are bounded by derivatives and terms that are not, so that the latter kind can be easily and reliably estimated from data, hence the ground-truth, while the first kind can be estimated by leveraging symbolic regression.\nSeveral synthetic datasets are employed for evaluation, many of which are simulated from equations that do not satisfy the linear combination assumption made by existing work.\n", "questions": "Please see my 'Weaknesses' section above.", "limitations": "The authors mentioned limitations. But in real-world scenario, there are other factors that might make applying this framework difficult. The first one is the sparsity of observations. The sensors are normally not well distributed and sometimes extremely sparse. So the estimate based on the zero-order information might not be reliable to start with. The second is the type of noise, which is normally unknown and needs to be estimated with the data together.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Strengths:\n\n1.\tA new framework for discovering PDEs, especially the ones with high-order derivatives and a lack of direct observations on the derivatives.\n\n2.\tA comprehensive evaluation on many different PDEs satisfying and beyond the assumption of PDEs made by previous work.\n\n3.\tThe comparison seems to show a better performance.\n\n4.\tGood exposition. The paper has a good balance between background and technical details.\n", "weaknesses": "I am in general in favour of this paper. However, it is not my area of expertise, so it would be good if they authors could clarify some questions here.\n\n1.\tReliance on symbolic regression.  I wonder to what extent the proposed framework has to rely on symbolic regression. This opens up several questions. (1) How inclusive or comprehensive does the dictionary has to be?  What if some key derivatives are not present in the dictionary?  (2) can the authors provide more details on the computational time in addition to E.2? It would be good to show the computation time in F.3, when the dictionary is gradually increased\n\n2.\tComparison. It is arguably true that derivatives are hard to measure in applications. But since the experiments are done using synthetic data, I wonder what if the observations of derivatives are available? Will the proposed method still outperform existing methods? I would assume the proposed framework still works to some extent, although not being able to make use of the observed derivatives.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer kC9e,\n\nWe appreciate your thorough feedback and kind remarks. Your questions have undoubtedly strengthened our work. In particular, we want to thank you for your inquiry about taking advantage of the derivative data should it be available. As D-CIPHER can in fact very naturally use such additional information, we have added a section in the appendix explaining how this can be achieved. We summarize the improvements we have made to the paper based on your review and we answer your questions below.\n\n### Actions taken\nWe have added the following sections to the appendix:\n1. Section \"How comprehensive does the dictionary $\\mathcal{Q}$ need to be?\" in Appendix F.\n2. Plot of the computation time as the dictionary is gradually increased in Appendix F.3\n3. Section \"How to take advantage of the observed derivatives?\" in Appendix F\n\n### Reliance on symbolic regression\nWe rely on the symbolic regression algorithm to find a closed-form derivative-free part. If the equation is suspected to have a derivative-free part equal to 0 or is linear in $u$ then the symbolic regression algorithm is not necessary and we can just add the term $u$ to the dictionary $\\mathcal{Q}$ (as a 0$^{\\text{th}}$ order derivative). This, however, only works in scenarios when the target equation is in a linear combination form. To discover anything more complicated, we need to rely on a symbolic regression algorithm.\n\n### How comprehensive does the dictionary $\\mathcal{Q}$ need to be?\nWe have attached a table (Table 1 in the attached PDF) that shows 17 different differential equations and what terms they need in a dictionary. All of them can be discovered with a dictionary of size of 10, $\\mathcal{Q} = (\\partial_t u, \\partial_x u, \\partial_t \\partial_x u, \\partial_t^2 u, \\partial_x^2 u, \\partial_t \\partial_x^2 u, \\partial_x^3 u, \\partial_x^4 u, \\partial_x(u^2), \\partial_x^2 (u^2))$, and 9 of those equations can be described with a dictionary containing just 4 terms, $\\mathcal{Q} = (\\partial_t u, \\partial_x u, \\partial_t^2 u, \\partial_x^2 u)$. It is thus likely that such small dictionaries are sufficient to discover most of the well-known equations. Many differential equations have similar derivative-bound parts and differ by derivative-free parts. Being able to discover *any closed-form derivative-free part* is what makes D-CIPHER stand out.\n\n### Computation time when the dictionary is gradually increased\nAs shown in Figure 1, CoLLie's computation time does not increase significantly when the dimensionality is increased. We have added in Appendix F.3 a plot that shows how the computation time increases when the dictionary is gradually increased. The plot can be seen in the attached PDF in Figure 1.\n\n### How to take advantage of the observed derivatives\nD-CIPHER can make use of the observed derivatives (if they are available) by adapting the dictionary. Consider a setting with a dictionary $\\mathcal{Q} = (\\partial_t u, \\partial_x u, \\partial_t \\partial_x u, \\partial_t^2 u, \\partial_x^2 u)$. If we happen to have the measurements of $\\partial_t u$ then we can introduce a new variable $v = \\partial_t u$ and change the dictionary to $\\mathcal{Q} = (v, \\partial_x u, \\partial_x v, \\partial_t v, \\partial_x^2 u)$. Note, we have performed experiments where the dictionary contains more than one dependent variable in Appendix F.2. With observed derivatives, we can also enlarge the space of Variational-Ready PDEs by allowing $g$ ($\\partial$-free part) to depend on $v$ as well.\n\nThank you once again for dedicating your time to reviewing our paper. We hope our responses have been satisfactory in addressing your queries. If any aspects still require additional explanation or if you have further questions, please let us know. We are eager to provide the necessary explanations."}, {"review_id": "iu8V6pMF7j", "reviewer": "Reviewer_Znhz", "summary": "This paper proposes a framework (D-CIPHER) to discover closed-form PDEs and ODEs. The framework is more general than some of the previously existing methods, and in particular can handle a class of PDEs defined as variational-ready PDEs in the paper. The empirical experiments evaluated the discovery performance on synthetic data for a set of different equations (showing both comparisons with existing methods on discovering linear combinations and results on discovering more challenging equations).", "questions": "See \"weaknesses\"", "limitations": "The very last paragraph discusses some potential limitations, but in my opinion it would be very helpful to have some \"negative examples\" in the paper, i.e. synthetic experiments where D-CIPHER fails and to explore the reasons of failure.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "Originality: From what I could tell, this paper proposes an original framework to discover broader classes of PDEs and ODEs. That being said, I'm not familiar with the literature in the area so I cannot fully speak to the originality.\n\nQuality: The empirical experiments, even though synthetic, appear to be thoughtfully designed and demonstrate notable improvements. \n\nClarity: This paper is very well-written with a good balance of technical details and general introductions. I enjoyed reading it even as someone outside of the ODE/PDE field.\n\n\n", "weaknesses": "\nSignificance: This is probably my biggest question for the paper. What types of real-world scenarios could D-CIPHER be applied to? The Discussion section briefly mentions finding heat and vibration sources and discovering population models and epidemiological models. However, at the level of the current discussion, these applications all sound very abstract. The paper would benefit from more grounding in concrete applications. If it's possible to add in experiments on real data, that would strengthen the paper. Even if not, the paper would still be improved with a more extensive discussion on how D-CIPHER could be applied in each of the relevant scenarios. For example, what would the form be (in step 1) based on prior knowledge? How would the fields be estimate (step 2)?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer Znhz,\n\nWe appreciate your thorough feedback and encouraging comments. We are particularly grateful for your suggestion on discussing how D-CIPHER can be applied in real-world scenarios. As a result, we have added two sections to Appendix F that discuss it in more detail. We firmly believe that these additions enhance the paper's overall quality and improve the accessibility of our approach for users.\n\n### Actions taken\nWe have added the following two sections to Appendix F:\n1. Impact on real-world problems\n2. D-CIPHER in practice.\n\n### Impact on real-world problems\nD-CIPHER is especially useful in discovering governing equations for systems with more than one independent variable. For instance, spatiotemporal data or temporal data structured by age or size. In particular, we envision D-CIPHER to be useful in modeling spatiotemporal physical systems, population models, and epidemiological models.\n\n**Spatiotemporal physical systems.** D-CIPHER may prove useful in discovering equations governing the oceans or the atmosphere. For instance, some places actively add or remove CO2 from the atmosphere. These “sources” and “sinks” are likely to be described by a $\\partial$-free part which D-CIPHER is specially equipped to discover. Similarly, with the ocean temperature where $\\partial$-free part can describe a heat source. Another area of application can be modeling seismic waves across the earth’s crust. Here the $\\partial$-free part can describe the vibration source (e.g., an earthquake).\n\n**Population models.** Population models can be used in agriculture to determine the harvest or for pest control to predict their impact on the crop. They have also been used in environmental conservation to model the population of endangered species. Population models have also been used in modeling the growth of cells to better understand tumor growth. Moreover, understanding the evolution of a population pyramid for a specific country may prove invaluable in ensuring its economic stability. As in all these scenarios, the rates of growth and mortality are likely to be described by $\\partial$-free part, D-CIPHER is uniquely positioned to discover such equations as an aid to human experts.\n\n**Epidemiological models.** Epidemiological models are crucial during a pandemic for better planning and interventions. For many diseases the rates of mortality and infection are age-dependent. Thus modeling the spread of disease using PDEs (rather than ODEs) might provide superior results. \n\n### D-CIPHER in practice\nWe have added a section to Appendix F discussing guidelines and things to consider while using D-CIPHER. In particular, it discusses the following points (a more detailed discussion is provided in the appendix).\n\n**The order of the differential equation.** One of the first considerations should be to choose the order of the differential equation $K$. For many dynamical systems, $K=2$ is sufficient unless we expect very complicated behavior. Then considering $K=3$ or even $K=4$ may be warranted. Note, that we show that D-CIPHER can discover a fourth-order PDE (Kuramoto-Sivashinsky equation) in Section 8.1.\n\n**Homogeneous equations.** Before searching through the whole space of closed-form $g$ (derivative-free parts), we can consider whether the equation we want to discover may be homogeneous. These experiments on the restricted search space can provide quick insights before searching through all closed-form derivative-free parts.\n\n**Terms in a dictionary.** For a given order of a differential equation $K$, it is a good idea to include all standard differential operators up to order $K$ acting on all the variables. For instance, for $K=2$ and $M=1+1$ we could choose $\\mathcal{Q} = (\\partial_t u, \\partial_x u, \\partial_t \\partial_x u, \\partial_t^2  u, \\partial_x^2 u)$. That allows to cover all linear PDEs with constant coefficients up to that order. To allow for non-linear PDEs we can include a term like $\\partial_t{u^2} = 2 u \\partial_t u$ that often describes advection (as in Burger's equation). \n\n**Dictionary steering when dealing with many dependent variables.** When we deal with a system of PDEs rather than a single PDE choosing a dictionary is increasingly important. As we explain in Appendix F.2, discovering whole systems of PDEs is very challenging and D-CIPHER is not designed to do so out of the box. However, we show how that can be done in certain situations. We can steer what kind of equations are discovered by choosing the terms in the dictionary.\n\n**Estimation algorithm.** Estimation algorithms make different assumptions on the data-generating process and should be chosen based on domain expertise. As we show in Appendix F.10, algorithms that produce smoother functions, such as Gaussian Process regression and cubic spline interpolation, tend to have good results. We can consider the advantages and disadvantages of these methods. For instance, Gaussian Process regression works very well for smooth signals. However, it is computationally intensive and might not perform well if the signal is not smooth enough (it has abrupt changes). Spline interpolation, on the other hand, is faster and more appropriate for less smooth signals, but it might introduce certain unwanted artifacts because of using cubic polynomials to interpolate the data.\n\n### Failure modes\nWe note that we not only discuss the *potential* limitations but actually show some of these limitations in our experiments. In particular, Figure 3 shows how D-CIPHER fails to discover the target equations if\n1. The noise becomes too high\n2. The sampling interval is too large\n3. The number of samples is too low\n\nThank you once again for dedicating your time to reviewing our paper. We hope our responses have been satisfactory in addressing your queries. If any aspects still require additional explanation or if you have further questions, please let us know. We are eager to provide the necessary explanations."}, {"review_id": "78e8d6TNU4", "reviewer": "Reviewer_ubhm", "summary": "The paper falls in the realm of data driven discovery of dynamical systems, PDEs to be specific. It proposes a framework for a class of PDEs which are termed as variation ready PDEs. These are claimed to be less restrictive than existing methods, which make stronger assumptions on the form of the PDE to be discovered and hence can't recover a significant population. A new optimization scheme, novel loss function and empirical evidence is provided to support the claims made in the paper.", "questions": "1) Can authors discuss the effect of the choice of test functions as B-splines? What about other options, how do they impact the results. A discussion on this ground would be helpful.\n2) Can authors shed light on the robustness of the proposed framework with respect to noisy and/or infrequent observations. These if theoretical can be taking into account sampling frequency, SNR, etc. \n3) How about parameterizing the test functions and making them learnable? Has this been tried, if not what do the authors feel in this regard. While test functions are fairly restrictive, there might still be a way to learn then or at the very least choosing them from a dictionary.", "limitations": "I don't see direct potential negative social impact. However, since this is a paper in the direction of ML for science and adjacent domains, it is quite possible that these methods when fully developed will have major impact. It will be good to have a word of caution regarding this, to make sure we acknowledge the vast amount of domain expertise already available in all scientific fields and not use ML as a tool to replace all human knowledge. \n\nAuthors have acknowledged technical limitations appropriately.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "In general the paper is well written, the scope of the paper is well thought. \n\nNotations are clear and introduced properly.\n\nIt address an important problem, the clear listing of challenges in the introduction is particularly impressive.\n\nThe variational loss function is novel", "weaknesses": "My major concern is regarding the way the problem is setup before the solution is proposed. Some of the terms introduced here although intuitive lack enough insight to make them more convincing. \n\nSection 7 is too short in the main paper to determine any novelty in the optimization scheme, this definitely needs to be presented better in the main paper, as this is claimed as a contribution.\n\nOne of the aspects mentioned in the paper is that of robustness from noisy or infrequent observations, I don't have see any evidence to support such claims.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer ubhm,\n\nWe appreciate your thorough feedback and encouraging comments. We summarize the improvements we have made to the paper based on your review and we answer your questions below.\n\n### Actions taken\n1. Elaborated on the requirements and other choices of testing functions in Appendix C.4.\n2. Provided more details of CoLLie in Section 7 and revised the contributions section.\n3. Added a discussion on the broader impact in Appendix F.\n\n### Problem setup\nWe are glad that you found the terms we define in Section 3 intuitive. We elaborate on their significance in Appendix F.7. \n\nTable 10, in particular, provides a concise summary of how the *evolution assumption* and the *linear combination* assumption impacts the optimization problem. As we explain in Section 3, we introduce *derivative-bound* and *derivative-free* parts of a PDE because they allow relaxing the two previously mentioned assumptions while still allowing for variational formulation (as not every PDE admits a variational formulation). This comes from a crucial observation that *no additional constraints* need to be put on the terms without derivatives to admit variational formulation.\n\nAs we write in Appendix F.7,  we believe the new notions complement the standard recognized PDE classes such as semi-linear, quasilinear, hyperbolic, etc. These standard classes were introduced predominantly to characterize the solving techniques or the properties of the solutions, whereas the notions we introduce relate to the difficulty of discovering such equations from data. Therefore we believe they offer substantial insight into this area.\n\n### CoLLie\nWe have added the following description of CoLLie to Section 7 and we modified the wording of the contribution section to indicate that CoLLie is a relatively minor contribution, compared with the main contribution of the paper, i.e., introducing new notions to describe PDEs from a discovery perspective, proposing a new general class of PDEs, a novel objective function, and a discovery algorithm.\n\n**Added to Section 7:**\n\nWe observe that this optimization problem is related to the one encountered in LASSO. Denote $\\mathbf{z}\\_0$ the solution that minimizes $||\\mathbf{A}\\mathbf{z}-\\mathbf{b}||\\_2^2$  (no constraints). If $||\\mathbf{z}\\_0|| > 1$ the problem is equivalent to finding $\\lambda$ (in the Lagrangian form of LASSO, Equation 48) such that the LASSO solution has the norm 1. Least Angle Regression (LARS) is a popular algorithm used to minimize the LASSO objective that computes complete solution paths. These paths show how the coefficients of the solution change as $\\lambda$ moves from $0$ to $\\lambda_{max}$ (from no constraints to effectively imposing the solution to be $\\mathbf{0}$). See Figure 6 in Appendix D2. CoLLie uses these solution paths to calculate the exact solution to the optimization problem. The case $0 < ||\\mathbf{z}_0|| < 1$ is harder as it corresponds to $\\lambda < 0$. CoLLie addresses this challenge by *extending* the solution paths generated by LARS beyond $\\lambda=0$ for $\\lambda<0$. We assume that the paths continue to be piecewise linear and keep their slope (Fig. 6, App. D2). CoLLie then uses these assumptions to efficiently find an approximate solution. We provide a detailed description of CoLLie in Appendix D.\n\n### Robustness to noisy and infrequent observations\nAs we show in Appendix F.5, estimating the derivatives is challenging and it becomes more difficult the higher the order of the derivative (see Figure 10). We thus show empirically that our method, which circumvents derivative estimation, performs better than the alternatives in scenarios with varying noise levels (Figure 2, Figure 3, Figure 4, Table 3, Table 8). We also show the impact of the frequency of observations and the number of samples in Figure 3. In all scenarios, our method is more robust than the alternative. We discuss the challenge of establishing theoretical error bounds in Appendix F.9.\n\n### Testing functions\nTesting functions should satisfy the following conditions:\n1. Be sufficiently smooth (at least $\\mathcal{C}^K$ for a $K^{th}$ order PDE)\n2. Compact support\n3. Derivatives can be computed analytically\n4. Orthonormal\n\nConditions 1 and 2 follow directly from Definition 3. Condition 3 is necessary because we do not want to estimate the derivatives of the testing functions. Condition 4 follows from the result obtained in [1] that suggests that these functions should be a subset of an orthonormal basis of L2 space.\n\nOther testing functions are possible and examining them constitutes an interesting research direction. In particular, piecewise polynomials as defined in [2] or various wavelets. Ideally, we would like to choose wavelets that form an orthonormal basis for the L2 space such as\n- Shannon wavelets - smooth ($\\mathcal{C}^{\\infty}$) but not compact\n- Meyer wavelets - smooth ($\\mathcal{C}^{\\infty}$) but not compact (better rate of decay than Shannon)\n- Daubechies wavelets - smooth ($\\mathcal{C}^{K}$ for a specified $K$) and compact but they do not have a closed-form expression.\n\nAnother interesting avenue of research would be to adapt the testing functions to the input data.\n\n### Impact\nWe have added a section that acknowledges that D-CIPHER is not designed to or capable of replacing human experts in scientific discovery and should be employed as a support tool as a part of a much broader scientific process.\n\n**References**\n1. Qian, Z., Kacprzyk, K. & van der Schaar, M. D-CODE: Discovering Closed-form ODEs from Observed Trajectories. ICLR (2022).\n2. Messenger, D. A. & Bortz, D. M. Weak SINDy for partial differential equations. Journal of Computational Physics (2021).\n\nThank you once again for dedicating your time to reviewing our paper. We hope our responses have been satisfactory in addressing your queries. If any aspects still require additional explanation or if you have further questions, please let us know. We are eager to provide the necessary explanations.\n\n\n\n\n"}]

Author rebuttal:
### Additional PDF

Table 1 in the attached PDF shows 17 different differential equations and what terms they need in a dictionary.

Figure 1 in the attached PDF shows how the computation time increases when the dictionary is gradually increased.

Meta-review: While the subject of the paper fell slightly out of the expertise of the reviewers, they collectively appreciated the strength of the contributions of the paper. The authors' response also addressed several raised concerns. Hence, I am recommending an acceptance.
Decision: Accept (poster)