Efficient Discovery of Dynamical Laws in Symbolic Form
Authors: Sören Becker
Michal Klein
Alexander Neitz
Giambattista Parascandolo
Niki Kilbertus
Venue: Submitted to ICLR 2023 (2023)
Keywords: Symbolic
ODE
Transformer

Abstract: We propose a transformer-based sequence-to-sequence model that recovers scalar ordinary differential equations (ODEs) in symbolic form from time-series data of a single observed solution trajectory of the ODE. Our method is efficiently scalable: after one-time pretraining on a large set of ODEs, we can infer the governing laws of a new observed solution in a few forward passes of the model. First, we generate and make available a large dataset of more than 3M ODEs together with more than 63M numerical solutions for different initial conditions that may serve as a useful benchmark for future work on machine learning for dynamical systems. Then we show that our model performs better or on par with existing methods in various test cases in terms of accurate symbolic recovery of the ODE, especially for more complex expressions. Reliably recovering the symbolic form of dynamical laws is important as it allows for further dissemination of the inferred dynamics as well as meaningful modifications for predictions under interventions.

Under review as a conference paper at ICLR 2023
EFFICIENT DISCOVERY OF DYNAMICAL LAWS IN SYMBOLIC FORM
An Anonymous authors Paper under double-blind review
ABSTRACT
We propose transformer-basee ordinary differential equations (ODEs) in symbolic form from time-series data of a single observed solution trajectory of the ODE, Our method efficiently scalable: after one-time pretraining on a large set of ODEs, we can infer the governing laws of a new observed solution in few forward passes of the model. First, we generate and make available a large dataset of more than 3M ODEs together with more than 63M numerical solutions for different initial conditions that may serve useful benchmark for future work on machine learning for dynamical systems. Then we show that our model performs better or on par with existing methods in various test cases in terms of accurate symbolic recovery of the ODE, especially for more complex expressions. Reliably recovering the symbolic form of dynamical laws is important it allows for further dissemination of the inferred dynamics as well meaningful modifications for predictions under interventions.
-oo-eueueee-eeeeeeee model that recovers scalar
INTRODUCTION
Science is commonly described as the "discovery of natural laws through experimentation and observation". Researchers in the natural sciences increasingly turn to machine learning (ML) to aid the discovery of natural laws from observational data alone, which often abundantly available, hoping to bypass expensive and cumbersome targeted experimentation. While there may be fundamental limitations what can be extracted from observations alone, recent successes of ML in the entire range of natural sciences provide ample reason for excitement. In this work, we focus on ordinary differential equations, a ubiquitous description of dynamical natural laws physics, chemistry, and systems biology. For a first order ODE = 0v/ac f t), call f (which uniquely defines the ODE) the underlying dynamical law. Informally, our goal is then to infer f in symbolic form given discrete time-series observations of single solution {yr yt,)) ) of the underlying ODE.
Contrary to "black-box-tedue box-techniques" such as Neural Ordinary Differential Equations (NODE) (Chen arguably opaque neural network, we focus
et al., 2018) that aim at inferring possible f specifically on symbolic regression. From the perspective of the sciences, a law of nature is useful insofar as it is more broadly applicable than to merely describe single observation. In particular, the reason to learn dynamical law in the first place is to dissect and understand it as well as to make predictions about situations that differ from the observed one. From this perspective, symbolic representation of the law (in our case the function f) has several advantages over block-box representations: they are compact and directly interpretable, they are amenable analytic analysis, they allow for meaningful changes and thus enable assessment of interventions and counterfactuals.
In this work, we develop Neural Symbolic Ordinary Differential Equation (NSODE), a sequence-to- sequence transformer to efficiently infer governing ODEs in symbolic form from single observed solution trajectory that makes use of massive pretraining. We first (randomly) generate total of >3M scalar, autonomous, non-linear, first-order ODEs, together with total of >63M numerical solutions from various (random) initial conditions. All solutions are carefully checked for convergence of the numerical integration. This dataset is unprecedented in both its scale and diversity and will be made publicly available alongside the code that was used to generate it.
We then devise NSODE, ie., numeric sequences of the form {(t., y: } directly to symbolic equations as strings, e.g., 22+1 64 cos (y) ", which is the prediction for f. This example directly highlights the benefit
peuueecc-to-sequeee
transformer that maps observed trajectories,

---

Under review as a conference paper at ICLR 2023
y f(y)
Operators
undry Brary
ODE solver
sample
thod numberr operatons
Model target
Linel 08
IEEE 7547
y+2. 3epxp coo
0
Symbols
2.3
sample
64
Model Inputs
Prior over ODEs
(g)
Data generation Training
Outputs (shifted right)
Two-hot encoding for numerical values
y+2 3texp( coa )
Transformer
Decodor
exp
-ommninn Cross-entropy loes
Next token prediction
Target
2.3
exp
Figure 1:
An overview illustration of the data generation (top) and training pipeline (bottom). Our dataset stores solutions in numerical (non- con-bbnariized) form on the entire regular solution time grid.
of symbolic representations in that the y² and cos(y) terms tell us something about the fundamental dynamics of the observed system; the constant 1. 64 will have semantic meaning in given context and we can, for example, make predictions about settings in which this constant takes different value. NSODE combines and innovates on technical advances regarding input representations and an efficiently optimizable loss formulation. Our model outperforms scalable existing methods in terms of skeleton recovery. especially on more complex expressions. While other methods still perform better on simple expressions found in (adjusted) existing benchmark datasets (which these methods have been tuned to) as well as a novel set of simple ODEs that we manually collected from different domains, they are typically orders of magnitude slower than NSODE
2 BACKGROUND AND RELATED WORK
Modeling dynamics and forecasting their behavior has a long history in machine learning. While NODEs (Chen et al. 2018) (with a large body of follow up work) are perhaps the most prominent approach, their inherent black-box character complicates scientific understanding
the observed
phenomena. Recent alternatives such tractable dendritic RNNs Brenner et al. 2022 and Neural Operators (Kovachki et al. 2021 et al. 202T set out to facilitate scientific discovery by combining accurate predictions with improved interpretability. A considerable advantage of these and similar approaches is their scalability to high dimensional systems as well as their relative robustness to noise, missing or irregularly sampled data (lakovlev et al. 2021 and challenging properties such as
multiscale dynamics [Vlachs et al. 2022 or chaos Park et al. 2022
Patel & OtL 2022).
Despite these advantages, we turn the focus different class of models in this paper and look at approaches that explicitly predict a mathematical expression in symbolic form that describes the observed dynamics. Such models are sense situated on the other side of the spectrum: while less scalable, the predicted symbolic expression is compact and readily interpretable so that dynamical properties can be analytically deduced and investigated. A recent benchmark study and great overview including deep learning- baaedd as well as symbolic models can be found in Gilpin 2021
Evolutionary algorithms. Classically, symbolic regression is approached through evolutionary algorithms such as genetic programming (Koza] 1993). Genetic programming randomly evolves
population of prospective mathematical expressions over many iterations and mimics natural selection by keeping only the best contenders across iterations, where superiority is measured by user-defined fitness functions (Schmidt & Lipson 2009). Process-bases modeling follows a similar approach but
includes domain sntwledge-informee constraints particular components of the system in order to reduce the search space to reasonable candidates [Todorovski & Dzeroski 1997 [Brdewwll et al.
2008 Simidjievski et al. 2020
2

---

Under review as a conference paper at ICLR 2023
Gradient-based search. More recently, symbolic regression has been approached with machine learning methods which exploit gradient information to optimize within the space of (finite) com- positions of pre-defined basis functions, [Bunton et al. 2016) builds on sparse linear regression to identify a linear combination of basis functions that yields the best fit for the observed data. This approach has inspired large body of follow-up work that generalize the idea partial observations
(Bakarji et al. 2022). parameterized functions (Lejarza & Baldea 2022). simultancous discovery of
coordinates Champion al. 2019 coordinate transformations that linearize the dynamics (Lusch et al. 2018 and partial differential equations (Rudy et al. 2017 These techniques often deploy regularizers and train one model for each set of observations. Once trained, the model itself represents the predicted symbolic expression, which can be read off the non-zero coefficients. This modeling principle is also employed by many other approaches that replace linear regression by neural networks with diverse sets of activation functions, both for differential equations
[ong et al. 2019 Liu et al. 2020 and non-differential algebraic equations (Sahoo et al. 2018).
Hybrid models. Alternatively, one can train model to directly output the symbolic expressions (for example as a string). Supervised learning with gradient-basee optimization for this approach is challenged by the formulation of differentiable loss that measures the fit between the predicted
symbolic expression and the observed data. Thus, prior work, mostly in the context of non-differential equation prediction, resorted reinforcement learning (Petersen et al. 2021 or combinations of
neural networks and evolutionary algorithms (Atkinson et al.] [201 9 Costa et al. 2021 1. A hybrid
approach combining gradient-free, frree, human intuition -guided heuristic search via genetic programming with neural network-based optimization has been presented for non-differential equations by Udrescu
al. 2020 and extended to dynamical systems by WWelbach et al. 202T
Sequence-to-sequemn models. The closest works to ours use pre-trained, attention- sequence-
bo-sequence models for symbolic regression of functional relationships (Biggio et al., 2021: Valipour
et al.
2021 Kamienny et al. 2022 Vastl et al.) (022) or (discrete) recurrence relations ( Ascoli
eet al. be both generated and evaluated on random inputs cheaply, resulting in essentially unlimited training data. Large data including ground-trun truth expressions in symbolic form allow for a differentiable cross-entropy loss based directly on the symbols of the expression, instead of the numerical proximity between evaluations of the predicted and true expression. While the cross-entropy loss works well for operators and symbols (e.g. +, exp, sin, x, y), a naive implementation is inefficient for numerical constants, e.g. 1. 452. As a workaround previous works resort to one of two strategies: 1) represent all constants with special <const> token when training the sequence-to-sseuennc model and predict only the presence of constant. Actual values are then inferred in a second, subsequent optimization step that minimizes Σi= Ily. der.ez... (I:) This second optimization procedure comes with substantial computational cost as constants have to be fit per inferred expression. In particular, we highlight that it does not transfer to inferring ODEs as would require to first solve the predicted ODE ý = f(y) to obtain predicted ( ) values that can be compared to the set of observations [ye While differentiable ODE solvers exist, optimizing constants this way is prohibitively expensive and typically highly unstable. 2) A second popular strategy consists in rounding constants within the range of interest so that they can be represented with finite number of tokens. This second strategy avoids a subsequent optimization step and enjoys clever encoding schemes with improved token efficiency yet represents values with inherent loss in precision. As an alternative, we develop representation based on a "two-hot" encoding which avoids subsequent optimization steps as well as rounding.
2022 They exploit the fact that symbolic expressions for (multi-variate) scalar functions can
We note that many previous symbolic regression methods have been described as "discovering natural laws". However, most of them learn fixed functional relationships from input-output pairs, whereas we seek to actually infer the underlying dynamic law that governed the behavior of the observed solution trajectory. A common obstacle to inferring dynamical laws is the lack of temporal derivatives in empirical measurements, which can hence not serve as regression targets. While derivatives can be estimated numerically, for example via finite differences, this requires dense, regularly sampled observations with a high signal-to-noise ratio to prevent fatal propagation of approximation errors. Alternative loss formulations that bypass unobserved time derivatives have been identified as an open challenge in symbolic regression for dynamical systems recently (Qian et al. 2022). With NSODE, we propose a solution that does not rely on temporal derivatives and thus avoids these complications.
3

---

Under review as a conference paper at ICLR 2023
3 METHOD
Problem setting. Given observations {(ts, k:) of trajectory y: [f1, [n] + R that is a solution of the ODE jj f(y), we aim to recover the function f in symbolic form-in our case as string. In this work, we focus on time-invariant (or autonomous) ODEs (i.e., f(y,t) are a good starting point for investigation as they are commonly studied and can be thought of as "evolving on their own" without external driving forces or controls, i.e. once an initial condition is fixed the absolute time does not directly influence the evolution. We explicitly assume that the observed system actually evolves according to an ODE in canonical form ÿ /= f(y) such that ʃ can be expressed in closed form using the mathematical operators seen during training (see Section 3.1]. To sum up. in this work we study the rich class of non-linear, scalar, first-order, autonomous ODEs, even though our approach is not fundamentally limited to this setting. In Appendix [A] we discuss extensions of NSODE to higher-order systems of coupled non-autonomous ODEs.
f()). Such settings
3.1 DATA GENERATION
Sampling symbolic expressions. To exploit large-scale supervised pretraining we randomly generate dataset of ~63 ODEs in symbolic form along with their numerical solutions for randomly sampled initial values. Since we assume ODEs to be in canonical form I f(). generating an ODE equivalent to generating a symbolic expression fly). We follow Lample & Charton (2019), who sample such an expression (1) as unary-binary tree, where each internal node corresponds to an operator and each leaf node corresponds to a constant or variable. The algorithm consists of two phases: (1) A unary-binary tree is sampled uniformly from the distribution of unary-binary trees with up to l e N internal nodes, which crucially does not overrepresent small trees corresponding to short expressions. Here the maximum number of internal nodes K is hyperparameter of the algorithm. (2) The sampled tree is "decorated", that is, each binary internal node is assigned a binary operator, each unary internal node is assigned unary operator, and each leaf is assigned variable or constant. Hence, we have to specify a distribution over the Nbin binary operators, one over the Nuna unary operators, a probability Psym E (0, 1) decide between symbols and constants, as well as distribution Pe over constants. For constants we distinguish explicitly between sampling an integer or a real value. Together with K, these choices uniquely determine distribution over equations f and are described in detail in Appendix B Figure depicts an overview of the data generation procedure.
The pre-order traversal of a sampled tree results in the symbolic expression for f in prefix notation. After conversion to the more common mathematical infix notation, we simplify each expression using the computer algebra system SymPy (Meurer et al. (017)) and filter out constant equations f(y) = as well as expressions that contain operators or symbols that were not part of the original distributior We call the structure modulo the value of the constants of such an expression a skeleton. Any skeleton containing at least one binary operator constant can be represented by different unary unary-binary trees. Vice versa many of the generated trees will be simplified to the same skeleton. To ensure diversity and to mitigate potential dataset bias towards particular expressions, we discard duplicates on the skeleton level. To further cheaply increase the variability of ODEs we sample Neoust unique sets of constants per skeleton. When sampling constants we take care not modify the canonical expression by adhering to the rules listed in Appendix B.I Our final dataset contains linear and non-lineat as well as homogeneous and inhomogeneous ODEs and we provide summary statistics about the distribution over equations in Appendix DJ Besides the number of internal nodes, we also report a simple yet common measure of complexity for each symbolic equation, which is the overall count of symbols (e.g., y. or constants) as well as operators in an expression.
Computing numerical solutions. We obtain numerical solutions for all ODEs via SciPy's 's interface (Virtanen et al. 2020) to the LSODA software package (Hindmarsh & Laboratory 1982) with both relative and absolute tolerances set to 10 LSODA consists of a collection of ODE solvers and implements strategy to automatically choose appropriate solver for the problem at hand (e.g.. recognizing stiff problems). We solve each equation on a fixed time interval / E [0, T] and store solutions on regular grid of N Ngrid points. For each ODE, we sample up to NA. initial values v() yo uniformly from (ymin yo maxx 2 While LSODA attempts to select an appropriate solver, numerical solutions still cannot be trusted in all cases. Therefore, we check the validity of solutions
With the exception of a unary -, which we do not discard.
Duee timcout per ODE, fewer solutions may remain if the solver fails for repeated initial value samples.
4

---

Under review as a conference paper at ICLR 2023
Table I: Overview of our model architecture.
Encoder
Decoder
architecture
BigBird'
BigBird
layers
6
6
heads
16
16
embed. dim.
512
512
forward dim.
2048
2048
activation
gelu
gelu
vocab, size
43
position enc.
Icarned
learned
parameters
23.3M
23.3M
We use full attention and chose BigBird (Zaheer al. 2020)) for its fast Huggingface implementation.
via the following quality control check: we use 9th order central finite differences to approximate the temporal derivative of the solution trajectory (on the same temporal grid as the proposed solution), denoted by yta, and filter out any solution for which U Vll e, where use 1.
3.2 MODEL DESIGN CHOICES
NSODE consists of an encoder-decoder transformer with architecture choices listed in Table 1 We provide visual overview in Figure
Representing input trajectories. A key difficulty feeding numerical solutions {u as input is that their range may differ greatly both within single solution as well across ODEs. For example, the linear ODE y = c for constant c is solved by an exponential y(t) = yo exp(ct) for initial value 0)) You which may span many orders of magnitude prevent numerical errors and vanishing or exploding gradients caused by the large range of values, assume each representable 64-bit float value is token and use its IEEE-754 encoding as the token representation (Biggio et al. 2021). We thus convert all pairs (1,, yt) to their IEFE-754 64 bit representations, channel them through a linear embedding layer, and then feed them the encoder.
a fixed time interval. To
Representing symbolic expressions. The target sequence (i.e., the string for the symbolic expression of is tokenized on the symbol-level. We distinguish two cases: (1) Operators and variables: for each operator and variable we include a unique token in the vocabulary. These tokens are one-hot encoded and passed through learnable embedding layer before their embedded representations are fed to the (2coder.(2) Numerical consiants: constants may come from both discrete (integers) as well as continuous distributions, as for example y* 2+1 64 + cos (y). Hence, unfeasible to include individual tokens "for each constant". Naively tokenizing on the digit level, i.e., representing real values literally as the sequence of characters (c.g., "1. 64 "), not only significantly expands the length of target sequences and thus the computational cost, but also requires a variable number of prediction steps for every single constant.
Instead, we take inspiration from Schrittwieser et al. (2200) and encode constants in two-hor fashion. We fix a finite homogeneous grid on the real numbers 21 < Ø2 < .... < Im for some m E N, which we add as tokens to the vocabulary. The grid range (æ1, and number of grid points m can be tuned for performance. Our choices are described in Appendix [B3] For any constant e in the target sequence we then find (l...., m 1} such that x; S Ii+1 and encode distribution supported on I;, Fi+1 with weights a, B such that ax; + 3z c. That is, the target in the cross-entropy loss for constant token is not strict one-hot encoding, but distribution supported on two (neighboring) vocabulary tokens resulting in a lossless encoding of continuous values in [21. ]. While this two-hot representation can be used directly in the cross-entee entsrpyy loss function and thus greatly facilitates training, it can not be passed directly through an embedding layer. For generic constant in the target sequence represented as a.r; + (Br++ 1. we thus instead provide the linear combination of the two embeddings a embed (x;) + embed (x111) ) decoder input.
Decoding constants. When decoding a predicted sequence, we check at each prediction step whether the arg max of the logits corresponds to one of the m constant tokens { Tm. ). If not, we proceed by conventional 1 1-hot decoding to obtain predicted operators and variables. instead the argmax corresponds to. for example. Fı also pick its largest-logit neighbor (x;-1 I suppose 1;+1), renormalize their probabilities by applying a softmax to all logits and use the resulting two probability estimates weights a, B. Constants are then ultimately decoded as ar, + 88441
5

---

Under review as a conference paper at ICLR 2023
3.3 EVALUATION AND METRICS
Sampling solutions. To infer a symbolic expression for the governing ODE of new observed solution trajectory {(4i, yi) H= all the typical policies such as greedy, sampling, or beam search are available. In our evaluation, we leverage computationally cheap forward passes to perform beam search with 1536 beams and report top-k results with & ranging from 1 to 1536. To pick from the candidate expressions provided by the beam search, one could for example only provide fraction of the trajectory as input to the model (say the first half) and then pick the candidate ODE whose numerical solution best predicts the trajectory on the remaining observation.
Metrics. We evaluate model performance both numerically and symbolically. For numerical evaluation we follow BBigii al. (2020): suppose the ground truth ODE is given by ÿ = [(v) with (numerical) solution y(t) and the predicted ODE is given by U f(y). To compute numerical accuracy we first evaluate if and ff on Neval points in the interval [min(((t)) max( (i.c., the interval traced out by the observed solution), which yields function evaluations gt = { Nayai and pred {6x Mara We then assess whether numpy- . a llllllls ose returns True well as whether the coefficient of determination R 2 0.999 Numerical evaluations capture how closely the predicted function approximates the ground truth function within the interval (min(y(t)), max (y(+))
However, a key motivation for symbolic regression is to uncover symbolic mathematical expression that governs the observations. Perhaps surprisingly, previous works on symbolic regression have paid little attention to directly testing for symbolic equivalence, which need not be implied by numerical fit. Testing for symbolic equivalence between ground truth expression (yy and a predicted expression fly) is unsuitable in the presence of real-valued constants as even minor deviations between true and predicted constants render the equivalence false. Instead, regard the predicted expression fly) to be symbolically correct if fyy and fly) can be made equivalent by modifying only the values of constants in the predicted expression fvv) This is implemented using SymPy's 's mat ch function. In order not to alter the structure of the predicted expression, constrain modifications of constants such that all constants remain non zero and retain their original sign. This definition is thus primarily concerned with the structure of an expression, rather than precise numerical agreement. Once the structure is known, the inference problem becomes conventional parameter estimation. We report percentages of samples in given test set that satisfies any individual metric (numerical and symbolic), as well as percentages satisfying symbolic and numerical metrics simultancously.
EXPERIMENTS
4.1 BENCHMARK DATASETS
We construct several test sets to evaluate model performance and generalization different settings.
testset-iv: Our first test set assesses generalization within initial values not seen during training. It consists of 5793 ODEs picked uniformly at random from our generated dataset but re-sampled initial values. We also employ the following constraints via rejection sampling: (a) All skeletons in testset-iv are unique. (b) As the number of unique skeletons increases with the number of operators, we allow at most 2000 examples per number of operators (with substantially fewer unique skeletons existing for few operators).
testset-constants: Our second test set assesses generalization within unseen initial values and constants. It consists of 2720 ODEs picked uniformly at random from our dataset (ensuring unique skeletons and at most 1000 examples per number of operators as above), but re-sampled intial values and constants.
testset-skeletons: In principle,
train NSODE on all possible expressions (using only the specified operators and number ranges) up to a specified number of operators. However, even with the millions examples in our dataset, we have by far not exhausted the huge space of possible
3 numpyy . al ase returns True if abs (a wise for elements a and b from the two input arrays. We use attolle--l 10 and rtol -0. 05; a corresponds to predictions, corresponds to ground truth.
b) (atol rtol abs (b) holds element-
*For observations ys and predictions ÿs we have R² 1 E, (.) ))

---

Under review as a conference paper at ICLR 2023
testset-iv
testset-constants
testset-skeletons
60
60
60
40
40
40
40
20
20
20
20
2"
2
2
2
2
2
2"
2
24
2
2
2
2"
2
2
2
2
2
2
2
2
2
2
2
topk
topk
topk
topk
(a) allelose
(b) R
2
0.999
(c) skeleton recovery
(d) skel. recov. & allelose
Figure 2: We evaluate numerical (allclose, R²) and symbolic (skeleton recovery) metrics well as combined skeleton recovery & allclose on testset-iv, testset-constanyy and testset-skeletons as k increases.
skeletons (especially for larger numbers of operators). Hence, our third test set contains 100 novel random ODEs with skeletons that were never seen during training.
testset-iv-166 This a subset of testset-iv motivated by the fact that most symbolic regression models we want to compare to require a separate optimization for every individual example, which was computationally infeasible for our testset-iv. For a fair comparison, we therefore subsampled up to 10 ODEs per complexity uniformly at random, yielding 163 examples.
To assess how NSODE performs on "real problems". we manually curated 12 scalar, non- non-linear ODEs from Wikipedia pages, physics textbooks, and lecture note from university courses on ODEs. These equations are listed in Table / in Appendix C We note that they are all extremely simple compared to the expressions in our generated dataset in that they are ultimately mostly low order polynomials, some of which with one fractional exponent.
testset-classic: To validate our approach on existing datasets we turn to benchmarks in the classic symbolic regression literature (inferring just the functional relationship between input-ouput pairs) and simply interpret functions as ODEs. In particular we include all scalar function listed in the overview in McDermott al. (2012) which includes equations from many different benchmarks Keijzer (2003); Koza ( [1993] 1994
;; [yy et al. (201T); Vladislavleva et al. (208). For example, we interpret the function f(y) = y² y² from Uy et al. (2011) 1) as an autonomous ODE
((t) fy((t) y(t)3 + y(t)² + y. which solve numerically for randomly sampled initial value as described before.
4.2 BASELINES
We compare our method to recent popular baselines from the literature (see Section 2 We briefly describe them including some limitations here and defer all details to Appendix L All baselines
work by explicitly fitting a separate regression function per individual ODE to map the observed samples of the solution trajectory to their temporal derivatives v() - i(), using the coefficient of
determination R² optimization objective. Since derivatives (t)) are typically not observed, approximate them via finite differences using the implementation available in PySindy (de Silva et al.) 2020). Hence, all these methods crucially rely regularly sampled and noise-free observations, whereas our approach can easily be extended to take those into account (see Appendix A).
Sindy (Brunton et al. 2016); Sindy builds a (sparse) linear combination of a fixed set of (non-linear) basis functions. The resulting Lasso regression is efficient however complex expressions such as nested or parameterized functions can only be approximated but not easily be represented in their precise symbolic form as all non-linear expressions have be explicitly added to the set of basis functions. We cross-validate Sindy over fairly extensive hyperparameter grid of 800 different combinations for each individual trajectory.
GPI (genetic programming): GPL(earn) maintains a population of programs each representing a mathematical expression. The programs are mutated for several generations to heuristically optimize a user defined fitness function. While not originally developed for ODEs, we can apply GPLeam on our datasets by leveraging the finite difference approximation. We use a population
gpl earn, rewdthedocs.ioo

---

Under review as a conference paper at ICLR 2023
Table 2: Comparing NSODE to the baselines on various benchmark datasets.
Dataset
Metric
NSODE
Sindy
GPLearn
AlFeynman
skel-recov
37.4
3.7
2.5
14.1
R² 2 0.999
24.5
31.9
3.7
49.7
iv-163
allclose
42.3
25.8
14.7
55.8
skel-recov & R² 0.999
15.3
3.1
1.8
13.5
skel-recov allciose
15.3
1.8
13.5
runtime [s]
5.4
0.4
29 +-22
1203.6
skel-recov
11.5
0
3.8
46.2
R² > 0.999
57.7
57.7
23.1
88.5
classic
allciose
80.8
57.7
30.8
88.5
skel-recov & R > 0.999
0
0
7.7
46.2
skel-recov allclose runtime
0
7.7
46.2
5.2
0.6
23 +22
1291.6
skel-recov
41.7
33.3
8.3
91.7
R 2 0.999
16.7
50
0.0
75
textbook
allclose
25
58.3
8.3
75
skel-recov & R 0.999
33.3
41.7
0
66.7
skel-recov & allclose
8.3
33.3
1.8
66.7
runtime [S]
23 +22
1267.
size of 1000 and report the best performance across all final programs. Compared to sindy, GPLearn is more expressive yet much slower to fit.
AlFennan (Udrescu & TTgmark 2020; Udrescu et al. 2020): AlFeynman is a physics-inspire approach symbolic regression that exploits the insight that many famous equations in natural sciences exhibit well-understood functional properties such as symmetries, compositionality, or smoothness. AlFeynman implements neural network based heuristic search that tests for such properties in order to identify a symbolic expression that fits the data. For every test sample AlFeynman computes a pareto front of solutions that trade off complexity versus accuracy. We report the best performance across all functions on the pareto front. Notably, AlFeynman performs quite an exhaustive search procedure such that running even on a single equation took on order of tens of minutes.
4.3 RESULTS
Model Performance. Figure shows NSODE's performance our testset-iv, testset-constants, and testset-skeletons according to our numerical and symbolic metrics as well as combined skeleton recovery and allclose sufficient test-time-coppute (i.e., considering many candidates) continuously improves performance. While capped k 1536 due to memory limitations, did not observe stagnation of the roughly logarithmic scaling of all performance metrics with k. This cannot be attributed to "exhaustion effects", where one may assume that all possible ODEs will eventually be among the candidates, because (a) the space of possible skeletons grows much faster than exponentially, and (b) the numerical metrics are extremely sensitive also to the predicted constant values in continuous domains. We verified this in an ablation experiment in Appendix
vary i in the top-k considered candidates of the beam search. Investing
As one may expect, performance decreases as we move from only new initial conditions, to also sampling new constants, to finally sampling entirely unseen skeletons. On testset-iv with 1: = 1536 we achieve about 50% skeleton recovery and still successfully recover more than a third skeletons of testset-skeletons with similar numbers for allclose. The fact that the combined metric (symbolic + numerical) is only about half of that indicates that numerical and symbolic are indeed two separate measures, none of which need imply the other. Hence, thorough evaluation of both is crucial to understand model performance symbolic regression tasks.
Comparison to baselines. In Table 2 we compare NSODE to all baselines (see Section 4.2 all our metrics (see Section (3.3) on all test sets (see Section|4] 4.1 using it 1536 in our beam search. We also include the average wallclock runtime per expression for each of the datasets. Since GPLeam
8

---

Under review as a conference paper at ICLR 2023
total
correct
10
10
10
count
count
count
count
20
20
20
20
complexity
complexity
complexity
complexity
(a) NSODE
(b) AlFeynman
(c) Sindy
(d) GPLcam
Figure 3: Correctly recovered skeletons by each method on basteet-iv-663 per complexity. AlFeynman and Sindy are mostly able recover some of the low complexity skeletons, while NSODE performs much better also on higher complexities. GPLearn fails to recover most skeletons.
often generates extremely long string expressions, it often takes SymPy up to half a minute to parse these expressions. We denote this extra time in gray and provide a more detailed version in Table g
First, we note that on our subsampled testset-iv-163, NSODE outperforms competing approaches in terms of skeleton recovery by a wide margin and also performs best in terms of joint skeleton recovery and numerical measures, which is a strong indication of actually having recovered the governing ODE accurately. By spending over 200x more time on its exhaustive heuristic search, A AlFeynman manages to outperform NSODE in terms of numerical accuracy (R² and allclose), which again does not necessarily correspond to correct identification of the ODE structure. Figure shows the number of skeletons recovered by each method given the complexity of equations, results for dataset-classic and datasct-textbook can be found in Appendix G While AlFeynman and Sindy recover some of the low complexity expressions, NSODE is the only method to also recover some of the more complex skeletons.
This also explains the strong performance on AlFeynman the overwhelmingly simple expressions in the testset-classic and testset-textbook. The heuristics deployed by AlFeynman are tuned to incorporate properties of known equations in the sciences. In particular, the algorithm explicitly attempts to fit a polynomial to the data which immediately corresponds to the correct ground truth expression in 8/ x 66.6% equations in testset-textbook and 12/26 46.2% equations testset- textbook. Together with the fact that it is solving somewhat simpler task (identifying functional relationships with manually provided derivatives from noise-free regularly sampled time-series), its strong performance on these particular test cases is to be expected. However, even these simple examples AlFeynman takes over 200x longer than our method, which in turn clearly outperforms Sindy and GPLearn in terms of skeleton recovery.
CONCLUSION
We have developed a flexible and scalable method infer ordinary differential equations (v) from a single observed solution trajectory. Our method follows the successful paradigm of large-scale pretraining of attention-based sequence-to-sequenee models on essentially unlimited amounts of simulated data, where the inputs are the observed solution {(t,, };) and the output is symbolic expression for f as string. Once trained, our method is orders magnitude more efficient than similarly expressive existing symbolic regression techniques that require a separate optimization for each instance and achieves strong performance terms of skeleton recovery especially for complex expressions on various benchmarks. Current limitations of our method mostly concern extensions to systems of high-order non-autonomous equations f. which we discuss in detail in Appendix A possible in principle with few to no modifications of NSODE, it is
While these extensions interesting direction for future to assess empirically how far this approach can be pushed. Besides the potential usefulness of our model in real-world scientific discovery and hypothesis generation, we also hope that our released dataset, code, and checkpoints will serve as a useful starting point for further research in this direction.
WWe note that the maximum complexity is not determined by K. maximum number of internal nodes in the tree alone, since also simplify all expressions using SymPy sometimes leading higher complexities.

---

Under review as a conference paper at ICLR 2023
6 ETHICS STATEMENT
This paper addresses fundamental research and as such its implications are difficult to forsee. However, all of our experiments are based non non-personalized data that we believe that datasets and model checkpoints can be publicly released without compromising private information.
REPRODUCIBILITY STATEMENT
To facilitate reproducibility we will release all datasets, source code our model and baseline implementations as well as our final model checkpoint. We also report all hyperparameter choices for our model in Section 3.2 and Appendix [B.3 as well as for the baseline methods Appendix E papre-parameters for data generation are described in Section 3.1 and concrete choices are listed in Appendix B.2
REFERENCES
Steven Atkinson, Waad Subber, Liping Wang. Genghis Khan, Philippe Hawi, and Roger Ghanem. Data-driven discovery of free-form governing differential equations. arXiv preprint 17X15:1910.05117, 2019.
Joseph Bakarji, Kathleen Champion, J Nathan Kutz, and Steven L Brunton. Discovering gov- erning equations from partial measurements with deep delay autoencoders. arXiv preprint arXiv:2201.05136, 2022.
Luca Biggio, Tommaso Bendinelli, Aurelien Lucchi, and Giambattista Parascandolo. A seq2seq approach symbolic regression. In Learning Meets Combinatorial Algorithms NeurIPS2020. 2020.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 936-945, PMLR, 18-24 Jul 2021. URL https:/prooceedngs..bh mlr. press/
v139/biggio2la.n la. at.
Manuel Brenner, Florian Hess, Jonas M Mikhaeil, Leonard F Bereska, Zahra Monfared, Po-Chen Kuo, and Daniel Durstewitz. Tractable dendritic mns for reconstructing nonlinear dynamical systems. In International Conference on Machine Learning. pp. 2292-222222222222222 2320. PMLR, 2022.
Will Bridewell, Pat Langley, Ljupčo Todorovski, and Sašo Džeroski. Inductive process modeling. Machine learning. 71(1):1-32, 2008.
Steven L. Brunton, Joshua L.. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedingy of the National Academy of Sciences, 1 113(15))3993223333, 2016. ISSN 0027-8424. doi: 10.1073/pnas. 15.1073/ 1517384113. URL
Kathleen Champion, Bethany Lusch, Nathan Kutz, and Steven L Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45): 22445-22451 2019.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 6572-6583, 2018.
Allan Costa, Rumen Dangovski, Owen Dugan, Samuel Kim, Pawan Goyal, Marin Soljačić, and Joseph Jacobson. Fast neural models for symbolic regression scale, 2021.
Stéphane D' Ascoli, Alexae-Ae Kamienny, Guillaume Lample, and Francois Charton. Deep symbolic regression for recurrence prediction. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 4520 4536. PMLR, 17-23 Jul 2022.
10

---

Under review as a conference paper at ICLR 2023
Brian de Silva, Kathleen Champion, Markus Quade, Dean-Christophe Loiseau, J Nathan 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): 2020.
William Gilpin. Chaos as an interpretable benchmark for forecasting and data-driven modelling. In Thirty-fifth Conference Neural Information Processing Systems Datasets and Benchmarks Track (Round 2), 2021. URI. httpps//opeeeevvi..ee//fooommnn/ ee Yit YJrf
A.C. Hindmarsh and Lawrence Livermore Laboratory. ODEPACK, Systematized Collection of ODE Solvers. Lawrence Livermore National Laboratory, 1982. URL https://booolsggggg...dbbbboo 9XWPmWEACAAI
wooglc.dc/boooko.c..... Pid=
Valerii lakovlev, Markus Heinonen, and Harri Lähdesmäki. Learning continuous time {pde)s from sparse data with graph neural networks. In International Conference on Learning Representations, 2021. URL http://peeneeiiiieeeetrrmommmdaa
Pierre-Alexandrs Alexandre Kamienny, Stéphane d' Ascoli, Guillaume Lample, and François Charton. End-to- end symbolic regression with transformers. arXiv preprint arXiv: 200., 2022.
Maarten Keijzer. Improving symbolic regression with interval arithmetic and linear scaling. In European Conference on Genetic Programming, pp. 70-82. Springer, 2003.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint arX4v12108.08881, 2021. doi:
John R. Koza, Genetic programming on the programming of computers by means of natural selection. Complex adaptive systems. MIT Press, 1993. ISBN 978-0-262-11170-6.
John Koza. Genetic programming II: automatic discovery of reusable programs. MIT press, 1994.
Guillaume Lample and François Charton. Deep learning for symbolic mathematics. In International Conference on Learning Representations, 2019.
Fernando Lejarza and Michael Baldea. Data-driven discovery of the governing equations of dynamical systems via moving horizon optimization. Scientific Reports, 12(1):1-15, 2022.
Zongyi Li, Miguel Liu-Schiaffini, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems, 2021.
Junyu Liu, Zichao Long, Ranran Wang, Jie Sun, and Bin Dong. Rode-net: learning ordinary differential equations with randomness from data, arXiv preprint arXiv:2006.02377 2020.
Zichao Long, Yiping Lu, and Bin Dong. Pde-net 2.0: Learning pdes from data with numeric- symbolic hybrid deep network. Journal of Computational Physics, 399:108992, 108925, 2019.
Bethany Lusch, Nathan Kutz, and Steven L Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature communications, 91)):1-10, 2018.
James McDermott, David R White, Sean Luke, Luca Manzoni, Mauro Castelli, Leonardo Vanneschi, Wojciech Jaskowski, Krzysztof Krawiec, Robin Harper, Kenneth De Jong, et al. Genetic program- ming needs better benchmarks. In Proceedings of the 14th annual conference on Genetic and evolutionary computation, pp. 791-798, 2012.
Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondiej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar. Sergiu Ivanov, Jason K. Moore. Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. Sympy: symbolic computing in python. PeerJ Computer Science, 3:e103, January 2017. ISSN 2376-5992. doi: 10.7717/peerj-cs. 103. URL https://b..oo///.777///pe......

---

Under review as a conference paper at ICLR 2023
Yonggi Park, Kelum Gajamannage, Dilhani I Jayathilake, and Erik M Bollt. Recurrent neural networks for dynamical systems; Applications to ordinary differential equations, collective motion, 2022.
and hydrological modeling. arXiv preprint
Dhruvit Patel and Edward Ott. Using machine learning to anticipate tipping points and extrapolate to post-tipping dynamics of non-stationary dynamical systems. arXiv preprint arXiv:2207.0052 2022.
Brenden K Petersen, Mikel Landajuela, T Nathan Mundhenk, Claudio P Santiago, Soo K. Kim, and Joanne T Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk- seeking policy gradients. In Proc. of the International Conference on Learning Representations, 2021.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela der Schaar. D-code: Discovering closed- odes from observed trajectories. In International Conference on Learning Representations, 2022.
Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science advances, ((4))=1602611 2017.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and
control. In International Conference on Machine Learning, pp. 4442 44500. PMLR, 2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324(5923):81-88) 2009.
Julian Schrittwieser, Ioannis Antonoglou, Thomas Hubert, Karen Simonyan, Laurent Sifre, Simon Schmitt, Arthur Guez, Edward Lockhart, Demis Hassabis, Thore Graepel, et al. Mastering atari, go, chess and shogi by planning with a learned model. Nature, 588(7839):604-609, 2020.
Nikola Simidjievski, L.jupčo Todorovski, Juš Kocijan, and Sašo Džeroski. Equation discovery for nonlinear system identification. IEEE Access, 6:29930-29943, 2020.
Steven H Strogatz. Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering. CRC press, 2018.
L.jupco Todorovski and Saso Dzeroski. Declarative bias in equation discovery. In Proceedings of the Fourteenth International Conference on Machine Learning, pp. 376-384, 1997.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: chysics-inspiredd method for symbolic regression. Science Advances, 6(16):eaay2631, 2020.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. arXiv preprint arXiv:2006. 10782, 2020.
Nguyen Quang Uy. Nguyen Xuan Hoai, Michael O'Neill, Robert McKay, and Edgar Galván- López. Semantically-biat crossover in genetic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Machines, 12(2):91-119, 2011.
Mojtaba Valipour, Maysum Panju, Bowen You, and Ali Ghodsi. SymbolicGPT: A Generative Transformer Model for Symbolic Regression. In Preprint Arxiv, 2021. URL httpss//aarriivorr//aa///
2106. 14131
Martin Vastl, Jonáš Kulhánek, Jifí Kubalík, Erik Derner, and Robert Babuška. Symformer: End-to- end symbolic regression using transformer-based architecture. arXiv preprint arX1v:2205.15764, 2022.
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, JJ Carey, ilhan Polat, Yu Feng. Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SeiPy 1.0: Fundamental Algorithms for Scientifie Computing in Python. Nature Methods, 17:261-272 2020. doi: 10. 092-019-
12

---

Under review as a conference paper at ICLR 2023
Pantelis R Vlachas, Georgios Arampatzis, Caroline Uhler, and Petros Koumoutsakos. Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelligence, (14),359-366, 2022.
Ekaterina J Vladislavleva, Guido F Smits, and Dick Den Hertog. Order of nonlinearity as a complexity measure for models generated by symbolic regression via pareto genetic programming. IEEE Transactions on Evolutionary Computationn 13(2):333-34, 2008.
Juliane Weilbach, Sebastian Gerwinn, Christian Weilbach, and Melih Kandemir. Inferring the structure of ordinary differential equations. arXiv preprint arXiv:2107.07:35, 2021.
Manzil Zaheer, Guru Guruganesh, Kumar Avinava Dubey, Joshua Ainslie, Chris Alberti, Santiago Ontanon, Philip Pham, Anirudh Ravula, Qifan Wang. Li Yang, et al. Big bird: Transformers for longer sequences. Advances in Neural Information Processing Systems, 33:17283- 17297, 2020.
A POSSIBLE EXTENSIONS
While we have focused exclusively on the huge class of explicit scalar, autonomous, first-order ODEs, we believe that our approach can scale also to non-autonomous, implicit, higher-order systems of ODEs.
Non-autonomous equations. Since our model is provided with time as inputs, it is capable of learning functions f(, t) depending on y and t explicitly. Hence, extending approach to non-autonomous ODEs is simply a matter of adding a symbol for t in our data generation.
Implicit equations. A generalization from purely explicit ODEs all ODEs of the form f(i. v) are possible but rely on the availability of appropriate ODE solvers for implicit equations. Sampling of this general form using the sampling algorithm described in Section B. straightforward as
we only need to include the additional variable ÿ in the set of candidate symbols and its associated distribution Psym- If we have an appropriate implicit numerical solver hand, we can supply the numerical solution trajectory to NSODE as before and use f for teacher forcing during training.
Irregular samples and noise. Due to the separate data generation and training phase, it is straight forward to train NSODE on corrupted input sequences, where we could for example add observation noise (need not be additive) to the {y: ) or randomly drop some of the observations on the regular observation grid to simulate irregularly sampled observations. We expect these models to still perform well and in particular much better than the baselines, since their targets are approximated derivatives, which are highly sensitive to noise and irregular samples.
Systems of equations. For a system of K first-order ODEs, the data generation requires two updates: first, we have to generate K equations (fi))e (I...... K). one for each component. In addition, each of those functions can depend not only on y (and t in the con-autonomous case), but K different components (x)( ke (1....KK This is casily achieved by allowing for K tokens y in the data generation. We can then simply augment the input sequence to the transformer to contain not only
(t,y) v) as for the scalar case, but (t, y(),, ( ((),....g{(k)) (K) for system of K variables. Finally, we have decide on a way to set the target for the transformer, i.e., how to represent the system of equations in symbolic form as a single sequence. A straight forward way to achieve this is to introduce special delimitation character, c.g.. to separate the different components. Note that the order in which to predict the K equations is dictated by the order in which they are stacked for the input sequence, hence this information is in principle available to the model. For example, in the case of a two-dimensional system, where y = (1) (2)) E R², we could have the model perform the following mapping
(1) (2) {l4, We Wi ) + "x1 (1) xy (2) x2(2) +cos (t) sin (y (1) + (2))) "",
where the separator delimits the two components of f(y, t) corresponding (() and 2(2)
Although the extension to systems is conceptually straightforward, we remark that the performance of the approach needs to be carefully evaluated in future work systems of ODEs show much larger variety of qualitative behaviors including chaos. One of the defining characteristics of chaos the sensitivity of solution trajectories to initial values (Strogatz) 2018), that is, even similar initial values
13

---

model with Under review as a conference paper at ICLR 2023
result in drastically different trajectories. Chaotic behavior may pose fundamental limitation for NSODE, which relies on supervised optimization on a representative training dataset, a challenging requirement in case of chaotic systems.
Higher-order equations. Finally, it is well known that any higher order (system of) ODEs can be reduced to first-order system of ODEs. Specifically, d-th order system of K variables can be reduced to an equivalent first-order system of d K variables. Hence, one can handle higher-order systems analogously as before with multiple separator tokens. One obstacle in this case is that when only observations of y() are given, one first needs to obtain observed derivatives to reduce higher-order system to a first-order system. These would turn have to be estimated from data, which suffers from the same challenges we have mentioned previously (instability under noise and irregular samples).
Finally, when we want to have a single model deal with higher-order equations of unknown order, or systems with differing numbers of variables, remains on open question how have the model automatically adjust to the potentially vastly differing input dimensions or how to automatically detect the order of an ODE.
B IMPLEMENTATION DETAILS
B. RULES TO RESAMPLE CONSTANTS
we generate ODEs as unary-binary trees, convert them into infix notation and parse them into canonical form using sympy. From each skeleton we then create up to 25 ODEs by sampling different values for the constants. When resampling constants we want to ensure that we do not accidentally modify the skeleton as this would additionally burden resolving potential ambiguities in the grammar of ODE expressions. Furthermore, we do not want to reintroduce duplicate samples the skeleton level after carefully filtering them out previously. We therefore introduce the following sampling rules for constants:
As described Section 3.1
3.1
1. Do not sample constants of value 0.
2. When the original constant in the skeleton is negative, sample negative constant, otherwise sample a positive constant.
3. Avoid base of in power operations 1 1.
4. Avoid exponent and -1 in power operations as x = x and x 1 1/2.
5. Avoid coefficients of value 1 and -1 as 1 x x and -1 -I
6. Avoid divisions by and 11 as x/ = x and x/ -x
B.2 DATA GENERATION
As discussed in the main text, the choices of the maximum number of internal nodes per tree K, the choice and distribution over Nai binary operators, the choice and distribution over Numa unary operators, the probability with which decorate a leaf with symbol Psym (versus constant with over constants uniquely determine the training distribution over ODEs f. These choices can be viewed as flexible and semantically interpretable tuning knobs to choose a prior over ODEs. For example, it may be known in a given context, that the system follows a "simple" law (small K) and does not contain exponential rates of change (do not include exp in the unary operators), and so on. The choice of the maximum number of operators per tree, how to sample the operators, and how to fill in the leaf nodes define the training distribution, providing us with flexible and semantically meaningful tuning knobs to choose prior over ODE systems for our model. We summarize our choices in Tables 3 to where '/ denotes the uniform distribution. Whenever a leaf node is decorated with a constant, the distribution over constants is determined by first determining whether to use an integer or a real value with equal probablity. In case an integer, sample it from Pint. and in case of a real-valued constant sample it from Preal shown in Table 3 Finally, when it comes to the numerical solutions of the sampled ODEs, we fixed the parameters in Table 6 for experiments.
1- Psym), ), and the distribution
14

---

Under review as a conference paper at ICLR 2023
We highlight that there is such thing "a natural distribution over equations" when it comes to ODEs. Hence, ad-hoc choices have to be made in one way or another. However, is important to note that neither our chosen range of integers the range of real values for constants are in any way restrictive as they can be achieved by appropriate rescaling. In particular, the model itself represents these constant values merely be non-numeric tokens and interpolates between those anchor tokens (our two-hot encoding) to represent continuous values. Hence, the model is entirely agnostic to the actual numerical range spanned by these fixed grid tokens, but the relative accuracy in recovering interpolated values will be constant and thus scale with the absolute chosen range. Therefore, scaling Piut and Preal by essentially any power of 10 does not affect our findings. Similarly, the chosen range of initial values (yo min yo non-restrictive as one could simply scale each observed trajectory to have its starting value lie within this range.
B.3 MODEL
For our Transformer model we choose the implementation of BigBird (Zaheer et al. 2020) available in HuggingFace. The model is trained 4 Nvidia A100 GPUs for 18 epochs after which
evaluate the best model based on the validation loss. We choose a batchsize of 600 samples and use linear learning rate warm-up over 10,000 optimization step after which we keep the learning rate constant at 10 For the fixed tokens that are used to decode constants, we choose an equidistant < F' 10 with m 21. This worked well empirically and using fewer or more tokens did not seem to improve model performance substantially. Finally, we use Newa 100 for the evaluation of our numerical metrics. We did not optimize architectural choices or hyperparameters for maximum performance and use the same choices in all experiments.
grid 10 21 < I2 <
While not relevant for our dataset as we check for convergence of the ODE solvers, we remark that the input-encoding via IEEE-754 binary representations also graciously represents special values such as nan or inf without causing errors. Those are thus valid inputs that may still provide useful training signal, e.g., "the solution of the ODE of interest goes to inf quickly".
C TEXTBOOK EQUATIONS DATASET
Table 7 list the equations we collected from wikipedia, textbooks and lecture notes together with the initial values that solved them for. We can also that almost all of these equations simplify to low-order polynomials.
D DATASET STATISTICS
We provide overview over the complexity distribution and the absolute frequency of all operators (after simplification) for all datasets in Figure We can see that our self-generated dataset covers by far the larges complexity whereas both complexities and operator diversity much lower for equations in the classic and textbook ODEs.
E BASELINES
We here describe more detail on the optimization of the baseline comparison models.
Sindy, We use the implementation available in PySindy (de Silva et al. 2020) and instantiate the basis functions with polynomials up degree 10 well all unary operators listed in Table 5 When fitting sindy to data we often encountered numerical issues especially when using high-degree polynomial the exponential function. To attenuate such issues we set the highest degree of the
Table 3: Parameter settings for the data generation.
parameter
K
Nbin
N Uma
Psym
Pint
Preal
value
5
0.5
u(1 10.
10) (0))
(0(())0,0))
15

---

Under review as a conference paper at ICLR 2023
Table 4: Binary operators with their relative sampling frequencies
operator
pow
probability
0.2
0.2
0.2
0.2
0.2
Table 5: Unary operators with their relative sampling frequencies.
operator
sin
cos
exp
log
probability
0.2
0.2
0.2
0.2
0.2
Table 6: Parameters for numerical solutions sampled ODEs.
parameter
Nconst
Niv
Ngrid
(vd 30
value
25
25
4
1024
(-5,5) 5)
Table 7: Equations of the textbook dataset.
Name
Equation f(x)
simplified
Yo
autonomous Riccati
0.6 0.6 1² + 0.1
0.6 2 2. 0.1
0.2
autonomous Stuart-Landau
-2.2/2 8¹ + 1.31
-1.1 y² 1.31
0.1
autonomous Bernoulli
-1.3 2.1 y222
-1.3 2.1 y²
0.6
compound interest Newton's 's law of cooling
0.1.y
0.1
-0.1 (y -
0.3 0.1
Logistic equation Logistic equation with harvesting
0.23 y (1 y)
0.23 (y y²)
0.23 y (1 - 0.33 3) 0.5
0.23 0.76 y² 0.5
Logistic equation with harvesting 2
2 y (1 v/3) 0.5
2. 0.66 y² 0.5
0.7
1/ (0.9 8 2.5). 1050,
Solow-Swan
7.2 y - 5.5
0.1
Tank draining
- (2/9)² v/y
00.11 y0.
Draining water
(0.53/4) V2- 9.81. (sin 1/ cos 1)2 y -1.
-0.67 /1..5
through funnel
velocity of body
99.81 - 0.9 9/8.2
-0.10 y 9.81
0.1
thrown vertically upwards
polynomials per sample to the highest degree present in the ground truth. Secondly, when numerical issues are caused by particular basis function, we discard this basis function for the current sample and restart the fitting process. We run a separate full grid search for every ODE over the following hyper-parameters and respective values (these all include the default values):
optimizer-threshold (np. 1 logspace (-5, 0, 10)): Minimum magnitude for a coefficient in the weight vector not be zeroed out.
optimizer-alpha ((0.001, 0.0025, 0.005, 0.01, 0.025, 0.05, 0.1, 0.2)): L2 regularizer on parameters.
finite differences order (2,,3, 5, 7, 9]): Order finite difference approximation.
maximum number of optimization iterations ([20, 100]): Maximum number of optimization steps.
For every ODE, sindy fit using solution trajectory the initial interval [o, 2] and validated on the interval (2, 4). The grid search thus results in ranking of models with different hyper-parameter configurations. Instead of evaluating only the performance of the best model, we report top-k performance across the ranked hyper parameter configurations. Sindy is computationally highly efficient yet complex expressions such as nested or parameterized functions can only be approximated (e.g. by polynomial basis) but not represented in their exact symbolic form.
16

---

Under review as a conference paper at ICLR 2023
training data
10a
count
106
count
10a
10
10
10c
5
10
15
20
25
30
add
cos
div
exp
log
mui
neg
Dow
sin
sqpt
siis
y
complexity
testset-iv
count
102
count
104
102
10
15
20
25
30
C
add
cos
div
exp
log
mu
mest
now
sin
sqru
mas
y
complexity
testset-const
10a
10a
count
10a
count
10
10a
10"
10-1
102
10
15
20
25
30
C
add
cos
div
exp
log
mul
neg
Dow
sin
sqri
suu
y
complexity
festset-skeleton
104
10a
count
10
count
10a
10
101
10-1
109
10
15
20
25
30
add
cOs
div
exp
log
mul
neg
Dow
sin
sqrt
suu
y
complexity
testset-iv-163
10
500
count
count
6
250
cos 1
0
5
10
15
20
25
30
add
div
exp
log
mull
neg
pow
sim
sgttt
complexity
tessset-classic
10
100
count
count
5
50
0
5
10
15
20
25
30
add
cos
exp
log
mal
neg
pow
sin
cart
sub
y
complexity
textsok-
5.0
40
count
count
2.5
20
0.0
add
div
mul
pow
sub
y
complexity
Figure 4: Distribution of complexity and operators for all datasets.
GPLearn. We instantiate GPLearn with constant range of (- 10) and all binary operators listed in Table and all unary operators listed Table 5 except for the exponential function which caused numerical issues. We keep the default papeeeparameters but grid search across the parsimony
17

---

Under review as a conference paper at ICLR 2023
coefficient (f0.0005, 0.001, 0.01, 0.05, 0.1, 0.2, 0.5, "auto")) which trades of fitness versus program length. We choose R² as fitness function.
AlFeynman. We use the AlFeynman implementation from http::/gttthbb.coo///o/lllA/ee// and run the algorithm with the following (default) paper--arametess:
Brute force try time: 60 seconds
Number of epochs for the training 500
Operator set: 14
Maximum degree of polynomial tried by the polynomial fit routine: 4
F ABLATION
The experimental results presented in Section 4.3 improves with increases in hyperparameter of the top-k sampling algorithm. While this result is intuitive, it also raises the question whether the model simply generates diverse set of random guesses that cover the full diversity of the testset. We investigated this question via the following ablation experiment: If we assume that the model only generates diverse set of random guesses, then we should be able to permute predictions across samples without loss in performance, c.g.
show that the performance measured by all metrics
(fi, (J2,
permute
(fi, fa.k))
(J2, U4,1, S4,x))
(J5r9a, Ufsro ... S6793,K)) (J5793, (f2, 12,
where fn the ground truth symbolic expression for sample and (f. fr, correspond to the generated predictions for sample n with top-k=K.
We implemented this experiment using the results obtained testset-iv for t 1536. Contrary to our assumption Table X shows that the resulting performance drops substantially in comparison to the original performance on all metrics. We can therefore conclude that our model is not simply generating diverse random guesses but actually learned a systematic mapping from input observations to symbolic expressions. Similarly, we can conclude that samples in the dataset are sufficiently distinct from each other, i.e. that the data generation process presented in Section 3.1 resulted diverse dataset.
Table 8: Comparison between original performance and performance after permuting predictions
across samples on testset-iv with k 1536. Results are given as percentages.
Metric
Original
Permuted
skel-recov
52.0
0.8
R² > 0.999
28.6
0.4
allclose
50.6
6.5
skel-recov & R² > 0.999
17.0
0
skel-recov & allclose
22.6
0.1
G DETAILED RESULTS
We provide comprehensive summary of performances of all models all datasets in Table 9
Figures 5 and further show again the number of correctly recovered skeletons by each method per complexity.
18

---

Under review as a conference paper at ICLR 2023
total
correct
NSODE
AlFeynman
7.5
7.5
count
5.0
count
5.0
2.5
2.5
0.0
0.0
15
20
10
25
30
complexity
complexity
Sindy
GPLearn
7.5
7.5
count
5.0
count
5.0
2.5
2.5
0.0
0.0
15
20
10
25
30
complexity
complexity
Figure Correctly recovered skeletons by each method on the classic benchmark dataset per complexity.
total
correct
NSODE
AlFeynman
4
count
count
2
10
15
20
25
30
10
15 complexity
20
25
30
complexity
Sindy
GPLearn
4
count
count
2
10
15
20
25
30
5
10
15
20
25
30
complexity
complexity
Figure Correctly recovered skeletons by each method on the textbook equation dataset per complexity.
19

---

Under review as a conference paper at ICLR 2023
Table 9: Detailed performance results for all methods on all (applicable) datasets.
Dataset
Metric
NSODE
Sindy
GPLearn AlFeynman
skel-recov
52.0
R² 0.999
28.6
lestset-iv
allclose
50.6
skel-recov & R 2 0.999
17
skel-recov & allclose
22.6
runtime [s]
5.3
skel-recov
45.6
R² 2 0.999
21.7
testset-constant
allclose
44.7
skel-recov & R² > 0.999
9.8
skel-recov & allclose
16.1
runtime [s]
5.3
skel-recov
19
R 0.999
12
testset-skel
allclose
33
skel-recov & R2 > 0.999
skel-recov & allclose
2
runtime [s]
5.3
skel-recov
37.4
3.7
2.5
R² 0.999
24.5
31.9
3.7
49.7
testset-iv-163
allclose
42.3
25.8
14.7
55.8
skel-recov & R² 0.999
15.3
3.1
1.8
13.5
skel-recov & allclose
15.3
3.1
1.8
13.5
runtime [s]
5.4
0.4
29
+22
1203.6
skel-recov
11.5
0
3.8
46.2
R > 0.999
57.7
57.7
23.1
88.5
classic
allclose
80.8
57.7
30.8
88.5
skel-recov & Rd > 0.999
0
0
7.7
46.2
skel-recov & allclose
7.7
46.2
runtime [s]
5.2
0.6
23 +22
1291.6
skel-recov
41.7
33.3
8.3
91.7
R² 0.999
16.7
50
0.0
75
textbook
allclose
25
58.3
8.3
75
skel-recov & R 2 0.999
33.3
41.7
0
66.7
skel-recov & allclose
8.3
33.3
1.8
66.7
runtime [s]
6
1
23 +22
1267.1
20


Reviews:
[{"review_id": "sBCV9kbls6", "reviewer": "Reviewer_UR3m", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this paper a transformer is trained to infer from a set of scalar time series observations the dynamical law behind these observations, in symbolic form. A large data set of equations and trajectories generated by these is generated, from which various test set scenarios are created. Performance is evaluated both in terms of the recovery of the “correct” symbolic expression, as well as numerically in terms of match between local time derivatives. The model is compared to various other approaches, and in general performs favorably, although it does not always outperform other models for reasons discussed.", "summary_of_the_review": "The authors may follow a useful, complementary direction, different from that most commonly engaged in scientific ML. However, I believe it is partly based on the wrong objectives if application in science is aimed for. It also falls far behind the methods based on universal approximation in expressivity, power, and scalability. Finally, Fig. 3 is not overly impressive (for any of the models tested), and the drop-off in success with model complexity hints at fundamental limits for this class of models regarding scalability to more realistic and higher-dim scenarios.", "strength_and_weaknesses": "Strengths:\n- yields a human-interpretable explicit symbolic form for the governing equation generating a given time series\n\nWeaknesses:\n- puts too much emphasis on an exact symbolic form, rather than on dynamically important properties\n- works only for very simple situations (scalar, noise-free, autonomous ...), unlike other methods developed in the dynamical systems field\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "hiJuUp8BB16", "reviewer": "Reviewer_LsD3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a transformer-based model to recover first order ordinary differential equations (ODE) in symbolic form, from points sampled from their trajectories at N regular time steps. In the paper, the trajectories are supposed to be noiseless, save perhaps for approximation error when solving the ODE, or rounding error of the values sampled on the trajectory. The ODE considered are of the form y' = f(y).\n\nThe model is a sequence to sequence 6 layer transformer (using HuggingFace implementation of BigBird as the base layer), with 512 dimensions and 16 heads (47M parameters overall). Training data are generated by sampling random functions f, using the \"tree algorithm\" introduced by Lample & Charton (2019), symplifying f with Sympy, solving the corresponding equation y' = f(y) with a numerical solver (LSODA), and sampling the corresponding trajectories for a number of initial values. The model loss is the cross-entropy of the sequence corresponding to f (after preorder enumeration of the corresponding tree). Beam search is used for generation.\n\nThe authors evaluate their model on three test sets: \n* solutions from the training set, but different initial values\n* solution skeletons from the training set, with different initial values, and constants\n* skeletons not found in the training set\nThey show that, with a large enough beam size, the model can reach significant accuracy (around 40%), both when numerically approximating f and recovering its skeleton. \n\nThey also provide (smaller) benchmark samples, and compare their model with three popular implementations : Sindy, GPLearn and AIFeynman. Their model outperforms Sindy and GPLearn, but AIFeynman, while over 200 times slower, is much more accurate. The authors suggest that this is due to differences between the training distribution and the benchmark samples (an out-of-distribution generalization problem).", "summary_of_the_review": "This is a very strong paper, on an important subject. It is clearly written, the experimental setting is very good, and the results are compelling. \n\nI recommend acceptance, and would be happy to increase my rating if the questions asked in the \"weaknesses\" sections are adressed during the rebuttal phase.", "strength_and_weaknesses": "**Strengths**\n\nThe paper adresses a hard and important problem of science: given numerical measurements over time, discover the underlying dynamics. \n\nThe paper is clearly written. The architecture and data generation techniques are clearly described and sound. \n\nThe experiments are well-thought, and support the claims. The comparisons with baselines are adequate, and I am very happy to see a negative result here (i.e. AIFeynman beating the model), together with a likely explanation of why this is so.\n\nOverall, a very solid paper, on an important subject. \n\n\n**Weaknesses** \n\nThere are no major weaknesses or flaws in the paper. I have a number of questions and suggestions for the authors.\n\n* **Limitations of the class of ODE**: at present, the only ODE considered have the form y'=f(y). A more general class of autonomous ODE would have the form f(y',y)=0, which could be generated with the same procedure, but will feature more difficult cases (i.e. singularities). The authors should mention this limitation, and possibly discuss it in appendix A. \n* **Simplifying training functions**: all the functions in the training and test set are simplified using Sympy. D'Ascoli et al. (https://proceedings.mlr.press/v162/d-ascoli22a.html, appendix C) observe that this is not necessary. On the other hand, it slows the generation procedure, and might make the model less robust since it reduce the diversity of the expressions used for training. Have the authors tried to train the model without simplification?\n* **Model ablation**: at present, the model uses a symmetric architecture, with 6 layers in the encoder and decoder. Some authors (e.g. https://arxiv.org/abs/2112.01898) have reported good results with asymmetric architectures. In this case, it seems that the decoder task is much harder, so maybe a shallow encoder and a deeper decoder would make sense. Kamienny uses a 16 layer decoder, and a 4 layer encoder. Have the authors experimented with such architectures? Also the models also use a fairly low dimension (512) and a large number of heads (16). Are there ablation studies justifying these choices? \n* **Scaling**: the models are relatively small (46 M parameters, 6 layers and 512 dimensions). On simpler problems, d'Ascoli uses 8 layers, and Kamienny 83M parameters. Experimenting with slightly larger models might improve performance, especially in the \"out of domain\" case of new skeletons (testset-skeleton).\n* **test-skeleton**: this test set is slightly smaller than the others. However, it provides the most practical assessment of model accuracy, since it measures model performance when the solution skeleton is not in the training set, a very likely case in practical situations. Could the authors provide a larger set when they make their dataset available? \n* **beam search**: experimental results suggest that very large beam sizes are needed for the model to perform correctly. D'Ascoli (https://proceedings.mlr.press/v162/d-ascoli22a.html, section 2.3) proposes a technique for ranking the hypotheses in the beam, and other techniques such as MCTS have been discussed as possible improvements over top-k perplexity beam search. Have the authors considered such approaches?\n* **benchmarks and out-of-domain generalization**: as the authors note, the better performance of AI-Feynman on benchmark sets might be explained by out-of-distribution generalization. These test sets feature much simpler functions that the ones seen at training. Have the authors considered retraining their model on a dataset of simpler functions (reducing the number of operators, and maybe changing the probabilities, in the algorithm they use), to provide for a better comparison (and support their interpretation of the benchmark results)?  ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "XuxJtYhglM", "reviewer": "Reviewer_qT9N", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose Neural Symbolic Ordinary Differential Equation (NSODE) a model mapping trajectories of numerical observations of the evolution of the behavior of a single variable through time to a symbolic expression that can generate those observations. They also present an approach to generating random symbolic expressions representing ODEs and simulations thereof starting from different initial conditions to train the model motivated by related work on symbolic regression for non-differential algebraic expressions.", "summary_of_the_review": "The presented approach while well written and even promising, it is severely limited and not correctly evaluated. Therefore, I lean towards rejection at this time.\n\n- Models only single dimensional ODEs.\nWhile the authors claim to discover dynamical laws, attempting to generalize to the general domain of science, the scope of model dynamical systems is severely limited. Physical systems tend to be more complex and models of such systems used in a broad range of scientific disciplines are represented as systems of differential equations.\n- Doesn't include approaches to inferring ODEs in the related work.\nIn the background and related work section the authors mostly discuss approaches to symbolic regression for non-differential algebraic expressions which address a different problem than the problem that the authors tackle with their approach. In the literature there exists a significant body of work addressing this problem. For example, Schmidt and Lipson 2009, or equation discovery for differential equations or process based modeling of dynamical systems. Furthermore, the only method mentioned in the related work claimed to address the problem of learning ODEs from data is NODE (Chan et al. 2018) which is in fact an approach to integrate (solve, simulate) differential equations.\n- Comparison to non-ODE baselines.\nRelated to the previous point, the authors compare to symbolic regression approaches that have not been designed for learning differential equations and yet interestingly outperform the proposed approach both with regards to several performance metrics and execution time.\n-Biased evaluation datasets.\nThe datasets that are generated using the approach can contain examples seen by the model and can be biased towards the presented approach. The main variability of examples in the dataset comes from the randomly selected initial conditions. Given the limitation to a single dimension, most of the equations are highly probable not to be sensitive to initial conditions and either converge to a single stable point including zero and infinity. Therefore multiple initial conditions can have very similar trajectories.", "strength_and_weaknesses": "Strengths:\n- Self contained. Presents a model and an approach to generate sufficient training data.\n- Comparative analysis.\n- Comparison on multiple datasets\n\nWeaknesses:\n- Models only single dimensional ODEs.\n- Doesn't include approaches to inferring ODEs in the related work.\n- Comparison to non-ODE baselines.\n- Biased evaluation datasets.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "2D5OvMnjrXB", "reviewer": "Reviewer_iM1D", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors present a Neural-Symbolic Ordinary Differential Equation (NSODE) that generates models straight out of the data. It is reported in the paper that once trained, the pipeline can produce models much faster than other techniques.", "summary_of_the_review": "The authors present a Neural-Symbolic Ordinary Differential Equation (NSODE) that generates analytical models that are simpler to handle. The presented experiments show that it is as good as current techniques.", "strength_and_weaknesses": "Strengths:\n1. The generation of symbolic models that use expressions commonly used by people to do analytical work is of tremendous importance.\n2. The proposed pipeline, once trained, is much faster than other techniques.\nWeaknesses:\n1. The authors state the Neural-Symbolic Ordinary Differential Equation (NSODE) is a symbolic model generator, as opposed to Neural Ordinary Differential Equations (NODE) (Chen et al., 2018) that generate opaque models. However, both develop symbolic models. The only difference is that the NSODE causes ones that use functions that we people understand and are helpful for analytical work. But, from a definition point of view, both approaches produce equally symbolic models.\n2. Even though the authors show that their \"model performs better or on par with existing methods,\" it still cannot resolve most of the presented problems, signaling that this work is inconclusive.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject