ODEFormer: Symbolic Regression of Dynamical Systems with Transformers
Authors: Stéphane d'Ascoli
Sören Becker
Philippe Schwaller
Alexander Mathis
Niki Kilbertus
Venue: ICLR 2024 spotlight (2024)
Keywords: symbolic regression
dynamical systems
differential equations
transformer
Primary area: general machine learning (i.e., none of the above)
TL;DR: We introduce ODEFormer, a transformer model able of inferring dynamical systems in symbolic form from observational data with state-of-the-art performance.

Abstract: We introduce ODEFormer, the first transformer able to infer multidimensional ordinary differential equation (ODE) systems in symbolic form from the observation of a single solution trajectory. We perform extensive evaluations on two datasets: (i) the existing ‘Strogatz’ dataset featuring two-dimensional systems; (ii) ODEBench, a collection of one- to four-dimensional systems that we carefully curated from the literature to provide a more holistic benchmark. ODEFormer consistently outperforms existing methods while displaying substantially improved robustness to noisy and irregularly sampled observations, as well as faster inference. We release our code, model and benchmark at https://github.com/sdascoli/odeformer.

Published as a conference paper at ICLR 2024
ODEFORMER: SYMBOLIC REGRESSION OF DYNAMICAL SYSTEMS WITH TRANSFORMERS
Stéphane d' d'Ascolii
Sören Becker*
EPFL
Helmholtz Munich
steephane.. dascol 1@gma 1 com
Munich Center for Machine Learning
TU Munich
soren. www.keeeeaaaaa.....o com
Alexander Mathis
Philippe Schwaller EPFL
Niki Kilbertus
EPFL
Helmholtz Munich
Munich Center for Machine Learning
TU Munich
ABSTRACT
We introduce ODEFormer, the first transformer able to infer multidimensional ordi- nary differential equation (ODE) systems in symbolic form from the observation single solution trajectory. We perform extensive evaluations two datasets: (i) the existing 'Strogatz' dataset featuring two-dimensional systems; (ii) ODEBench, collection of one- to four-dimensional systems that carefully curated from the literature to provide a more holistic benchmark. ODEFormer consistently outperforms existing methods while displaying substantially improved robustness to noisy and irregularly sampled observations, as well as faster inference. We release our code, model and benchmark at https://github.co http:///hhubooom//ssc////oooooo/
INTRODUCTION
Recent triumphs of machine learning (ML) spark growing enthusiasm for accelerating scientific discovery (Davies al., 2021: Jumper al., 2021: Degrave et al., 2022). In particular, inferring dynamical laws governing observational data is an extremely challenging task that is anticipated to benefit substantially from modern MI. methods. Modeling dynamical systems for forecasting, control, and system identification has been studied by various communities within MI.. Successful modern approaches are primarily based on advances in deep learning, such neural ordinary differential equation (NODE) (see Chen al. (2018) and many extensions thereof). However, these models typically lack interpretability due to their black-box nature, which has inspired extensive research explainable ML of overparameterized models (Tjoa & Guan, 2020; Dwivedi et al., 2023). An alternative approach is to directly infer human- readable representations of dynamical laws.
This is the main goal of symbolic regression (SR) techniques, which make predictions in the form of explicit symbolic mathematical expressions directly from observations. Recent advances in SR make it promising alternative to infer natural laws from data and have catalyzed initial successes in accelerating scientific discoveries (Aréchiga et al., 2021; Udrescu & Tegmark, 2021; Butter et al., 2021). So far, SR has most commonly been used to infer function g(x) from paired observations (x, (x)) we call this functional SR. However, in many fields of science, understanding a system involves deciphering its dynamics, i.e., inferring function fx) governing its evolution via ODE i f(x)- we call this setting dynamical SR. The task is then to uncover f from an observed solution trajectory (1, (()), where observations of (t) may noisy and times /. may be irregularly sampled.
Contributions. In this work, we introduce ODEFormer, the first Transformer trained to infer dynamical laws in the form of multidimensional ODEs from observational (noisy, irregularly sampled) data. It relies large-scale training of a requence-to-sequee transformer diverse synthetic data, leading to efficient and scalable inference for unseen trajectories. Faced with the lack of benchmarks
"qqual contribution.

---

Published as a conference paper at ICLR 2024
Lang unge death, mnddd language
Landau Oupical Lime scals tau 1]
Lotka-Nalhrra simple Wikipodia)
non-lnnar cscillation)
Moxwwell Bloch equations dynamics)
Loronz equattons in perooiic regtme
Btnocular stralry model with adaptatton loscLiations)
SEIR infaction model (proportions)
0
10
10
0
5 10
5 10
Figure I: ODEFormer predicts dynamics from noisy and irregularly sampled trajectories in 1-4 dimensional systems. Ground truth trajectories (thick lines) are corrupted (5% noise) and unevenly sampled (50% of the equally spaced points are dropped uniformly at random), yielding the observed data (dots) from which ODËFormer infers a symbolic ODE. Integrating the predicted ODE from the original initial values gives trajectory (thin lines) that closely approximates the observed data.
for dynamical SR (the only existing one, called "Strogatz dataset" (La Cava al., 2021), contains only seven two-dimensional systems, and is not integrated with sufficient precision (Omejc et al., 2023)), we also introduce ODEBench, a more extensive dataset of 63 ODEs curated from the literature, which model real-world phenomena in dimensions one to four (see Appendix A for details). On both benchmarks, ODEFormer achieves higher accuracy than existing methods for irregular, noisy observations while being faster than most competitors.
Problem setting and overview. We assume observations {(ta, 24) [E(j. of solution x(t) of = f(x), for some f D + RD, where x; can be noisy observations of xr(4) with irregularly sampled f.. Here, task is to infer f in symbolic form from {(t, x)) hEgNN. As illustrated in Figure 2, ODEFormer is based on large-scale pre-training on "labelled" examples that consist of randomly generated symbolic mathematical expressions as the prediction target ʃ and a discrete solution trajectory {(4, (i),)))) as input, obtained by integrating i = fzz) for a random initial condition.
use the notation [N] == {1, (1.....4 NN and 4 for the temporal derivative. The
2 RELATED WORK
Modeling dynamical systems. We briefly mention some relevant cornerstones from the long history of modeling dynamical systems from data, each of which has inspired a large body of follow-up work. Neural ODEs (NODE) (Chen et al., 2018) parameterize the ODE f by neural network and train it by backpropagating through ODE solvers (either directly or via adjoint sensitivity methods). NODEs require domain knowledge but only represent the dynamics as generally overparameterized black-box model void of interpretabiliity. Assuming prior knowledge about the observed data, physics- informed neural networks (Raissi et al., 2019; Karniadakis al., 2021) aim model dynamical systems using neural networks regularized to satisfy a set of physical constraints. This approach was recently extended uncover unknown terms in differential equations (Podina et al., 2023). In this work, we aim to infer interpretable dynamical equations when no domain knowledge is available.
Approaches to symbolic regression. While it is possible to generate complex analytical expressions that fit the observations, an unnecessarily lengthy function is often impractical. Symbolic regression seeks balance between fidelity to the data and simplicity of form. Therefore, predictions are typically also evaluated in terms of some "complexity" metric. Because the symbolic output makes difficult to formulate differentiable losses, SR has traditionally benefitted comparably little from advances in autodifferentiation and gradient-based optimization frameworks. The dominant approach has thus been based on evolutionary algorithms such as genetic programming (GP) (La Cava et al., 2016b; 2018; Kommenda et al., 2020; Virgolin et al., 2021; Tohme et al., 2022; Cranmer, 2023), optionally guided by neural networks (Mundhenk et al., 2021; Udrescu & Tegmark, 2020; Costa et al., 2021) and recently also employing reinforcement learning (Petersen et al., 2020) see (La Cava et al., 2021; Makke & Chawla, 2022) for reviews. Most of these approaches require a separate optimization for each new observed system, severely limiting scalability.
2

---

Published as a conference paper at ICLR 2024
(+, 9783, -0, (+, 1735, E-31,
(c, 8134, E-31, 1+, 5992, E-9).
FFN- -
Encodier
Decoder
[+, 1034, E-4), (+, 4305, E-4),..
Gunb
time
Target $20.
Observations
Tokenization
Embedding
Transformer
Cross-Entropy loss
Figure 2: Model overview. ODEFormer maps (noisy) observations to symbolic sequences that represent (the right-hand side of) an ordinary differential equation (ODE). Model input includes representations of all state variables (blue & orange) and time (gray). The model is optimized via cross-entropy loss between the softmax distribution over the model vocabulary (gray) vs. the one-hot encoded true target symbol (green). Target sequences are predicted in prefix notation. The illustration includes the first three predictions [+, Z2, 21] for the target sequence [+, 21, xz2._ This example showcases that the loss may be high despite mathematical equivalence between predicted and target sequence this does not impact training practice.
Transformers for symbolic regression. With the advent of transformer models (Vaswani et al., 2017), efficient learning of sequence-to-sequeeee tasks for broad variety of modalities has become feasible. Paired with large-scale pre-training on synthetic data, transformers have been used for symbolic tasks such as integration (Lample & Charton, 2019), formal logic (Hahn et al., 2021), and theorem proving (Polu & Sutskever, 2020). Few recent works applied them to functional SR (Biggio et al., 2021; Valipour et al., 2021; Kamienny et al., 2022; Vastl et al., 2022) obtaining comparable results to GP methods, with key advantage: after one-time pre-training, inference is often orders of magnitude faster since no training is needed for previously unseen systems. Landajuela et al. (2022) recently proposed a hybrid system combining and leveraging the advantages of most previous approaches for state of the art performance on functional SR tasks.
Symbolic regression for dynamical systems. In principle, dynamical SR, inferring f from (t, (t)), can be framed as functional SR for ((f)),(()) pairs. However, when transitioning from functional to dynamical symbolic regression, key challenge is the absence of regression targets #(t) since temporal derivatives are usually not observed directly. A common remedy is to use numerical approximations of the missing derivatives as surrogate targets instead. This approach has been employed in the GP community (Gaucel et al., 2014; La Cava et al., 2016a; Kronberger et al., 2020; Quade al., 2016) and also key to the widely used SINDy (Brunton al., 2016) algorithm which performs sparse linear regression on manually pre-defined set of basis functions. While SINDy is computationally efficient, its modeling capacity is limited linear combinations of its basis functions. Like in functional SR, neural networks have also been combined with GP for dynamical SR (Atkinson et al., 2019), and the andie--and-conque strategy by Udrescu & Tegmark (2020) has also been extended to inferring dynamical systems (Weilbach et al., 2021). Omeje et al. (2023) recently introduced ProGED, which performs dynamical SR via random search of candidate equations, constrained by probabilistic context-free frree grammars (PCFGs). However, this approach requires knowledge about the ground truth as parameters of the PCFGs need to be carefully tailored to each problem. Finally, Qian et al. (2021) recently presented innovative loss function to mitigate finite difference approximations which are error-prone in the presence of signal corruptions. In the realm of transformers, d' Ascoli al. (2022) infer one-dimensional recurrence relations in sequences of numbers, which as discrete maps are closely related to differential equations. Most related to our work. Becker et al. (2023) explore a transformer--asse approach dynamical SR for ODEs (for a detailed comparison, see Appendix H). However, their method is limited to univariate ODEs. Such systems exhibit extremely limited behavior, where solution trajectories can only either monotonically diverge or monotonically approach a fixed value not even inflections let alone oscillations are possible (Strogatz, 2000). In this work, we tackle the important yet unsolved task of efficiently discovering arbitrary non-linear ODE systems in symbolic form directly from data in multiple dimensions, without assuming prior knowledge of the ground truth.
Theoretical identifiability. Traditionally, the identifiability of ODEs from data has been studied in case-by-case analysis for pre-defined parametric functions (Áström & Eykhoff, 1971; Miao et al., 2011; Villaverde et al., 2016; Hamelin et al., 2020). Recent work made progress on identifiability within function classes such as linear (in parameters) autonomous ODE systems (Stanhope et al., 2014; Duan et al., 2020; Qiu et al., 2022) or, in the scalar case, even non non-parametric classes such as analytic. algebraic, continuous, or smooth functions (Scholl et al., 2023). As it stands, it is difficult to gain practical insights from these results for our setting as they do not conclude whether non- non-linear ODEs can practically be uniquely inferred from data Appendix B for details.
3

---

Published as a conference paper at ICLR 2024
3 DATA GENERATION
Our method builds on pretraining on a large dataset of ODEs which is assembled as follows.
Generating ODEs. For D-dimensional ODE f. we independently sample the D component
functions fi....: fo as random unary-binary trees following the method of (Lample & Charton, 2019), where internal nodes represent mathematical operators and leaves represent constants or variables. In our specific procedure, we first sample the system dimension D uniformly from [Dmax] for fixed max E N and then perform the following steps for each component function:
1. Sample the number of binary operators b uniformly from [bmax]] for fixed bmax N.
2. Sample a binary tree with non-leaf nodes, following the procedure of Lample & Charton (2019).
3. Decorate each non-leaf node with a binary operator sampled from P(+) = 3/4 and P(x) 1/a.
4. For each leaf in the tree, sample one of the variables I; for E [D]-
5. Sample the number of unary operators uniformly from [Mmaax]] for fixed Vmnax E N.
6. Iteratively with u repetitions, select a node whose subtree has a depth smaller than 62 and insert a new node directly above. Populate the new node with unary operator that is sampled uniformly at random from {x + sin(x), . + x== x + x² )²).
Convert the tree into a mathematical expression via preorder traversal (Lample & Charton, 2019). 8. Finally, prepend coefficient to each term and wrap the argument of any unary operator in an affine transformation x + a x + b. Coefficients and constants of affine transformations are sampled independently from a log-uniform distribution on [Cmin> Cmax] with Cmin, Cmax R.
Due to random continuous constants (and initial conditions), we almost surely never sample function twice. In our experiments, we use D Dmax= 6, bmax
5. V/max 3, (Cmin: Cmax) = (0.05, (0).
Integrating ODEs. Once the function f is generated, we sample an initial condition zo (0, 7LD) for a fixed y E R= and identity matrix II, and integrate the ODE from
= T using the numerical solver scipy. int integrate. solve_ ivp provided by SciPy (Virtanen et al., 2020). function defaults to an adaptive 5th order explicit Runge Kutta method with 4th order error control (Dormand & Prince, 1980) as well as relative and absolute tolerances of 10-3 and 10-6 respectively. To avoid the sampling grid to be informative of signal variation, we use the solver's internal option to interpolate the adaptive step-size solution to obtain the final solution cquidistant grid of N sampling points, where N ~ 550, 2003. When integration fails, i.c., when the solver throws error, returns unsuccessfully, takes longer than one second, we simply discard the current sample. Section 4 explains in detail that we can fix T and y during training without loss of generality, we can rescale observations inference time. Hence, we use 7 1 and T 10.
The solve_ive
Filtering data. Under the distribution over functions / defined implicitly by our generator, a substantial fraction of sampled ODEs (and initial conditions) leads to solutions where at least one component diverges over time. Continued divergence over long time spans is typically deemed "unphysical". Among the trajectories that remain finite, we observe that again a substantial fraction swiftly converges towards fixed point. Although these solutions may be realistic, their dominance hampers diversity in our dataset. Also, stable constant observations over long times spans are arguably rarely of interest. Hence, we use the following heuristics to increase diversity of the generated ODEs:
If any variable of the solution trajectory exceeds fixed threshold (10² in our experiments), we discard the example. This amounts to filtering out divergent systems.
If the oscillation of all component functions over the last quarter of the integration range is below certain threshold (10-3 in our experiments), we discard the example with a probability of 90%. This filters out a majority of rapidly converging systems.
MODEL, TRAINING, AND INFERENCE
ODEFormer is an encoder-decoder transformer (Vaswani et al., 2017) for end-to-end dynamical SR as illustrated in Figure 2. The model comprises 16 attention heads and an embedding dimension
Subtractions and divisions are included via multiplication with negative numbers and the unary operator + respectively. It has been argued that divisions appear less frequently than additions and multiplications "typical" expressions (Guimerà al., 2020).
"This aims at avoiding deeply nested uninterpretable expressions, which often occur in GP-based SR.
"The oscillation of function h (a, b] R is given by supac(a.M h(x) -nn h(z).
4

---

Published as a conference paper at ICLR 2024
of 512, leading to a total parameter count of 86M. As observed by Charton (2022), we find that optimal performance is achieved in an asymmetric architecture, using 4 layers in the encoder and 16 in the decoder. As measurement time is explicitly included in the imputs, we remove positional emheddings from the encoder. ODFFormer is trained on the data set of 50M samples generated as described in Section 3 following established principles, with details given in Appendix C. Whereas the transformer itself follows the standard architecture, we adapted tokenization and embedding to the particular task, i.e., inference of symbolic functions from multivariate observations.
Tokenizing numbers. Since numeric input trajectories as well as symbolic target sequences may contain floating point values, we need an efficient encoding scheme that allows the infinite number of floats to be sufficiently well conserved by fixed-size vocabulary. Following Charton (2022), each number rounded to four significant digits and disassembled into three components: sign, mantissa and exponent, each of which is represented by its own token. This tokenization scheme condenses the vocabulary size required to represent floating point values to just 10203 tokens (+, 0, .... 9999, E-100 1000,-. .... E)C0) and works well in practice despite the inherent loss of precision. We also experimented with three alternative representations: (i) two-token encoding, where the sign and mantissa all merged together, (iii) a two-hot encoding inspired by Schrittwieser et al. (2020) and used by Becker et al. (2023), which interpolates linearly between fixed, pre-set values to represent continuous values. These representations have the advantage of decreasing sequence length, and (iii) has the added benefit of increased numerical precision for the inputs. Since all three alternatives led to worse overall performance, we used the three token representation (sign, mantissa, exponent).
merged together, (ii) one-token encoding where sign, mantissa and exponent
Embedding observed trajectories. For D dimensional ODE system, the above tokenization scheme represents data of a single observed time- tine-point (t;,,;) using R(++1) x 3 tokens. Each such token is fed separately to an emhedding layer and the resulting token embeddings are concatenated to give a combined embedding vector in [(((41)x3)xden for each observed time-point i. We want single ODEFormer model to flexibly handle systems of different dimensionalities (up to D Duax= = 6). In case the observed system has dimension D < Dmaa we hence fill vacant dimensions via seero-padding Following Kamienny et al. (2022), the resulting embedding vector subsequently processed by a 2-layer fully-connected feed-forward network (FFN) with Sigmoid-weighted linear unit (SiLU) activation functions (Elfwing al., 2018), which reduces the dimension of the final representation of single time point to demb dimensions. This embedding process allows a single trained model to flexibly handle input trajectories of variable lengths and ODE system dimensions.
Encoding symbolic functions. In addition to the tokens used to represent floating point values as described above, the vocabulary of the decoder includes specific tokens for all operators and variables. Importantly, the decoder is trained on expressions in prefix notations to relieve the model from predicting parentheses (Lample & Charton, 2019). For D-dimensional systems, we simply concatenate the encodings of the D component functions, separated by special token "I". With these choices, the target sequence for an exemplary 2-dimensional ODE given by component functions fi cos 21) + I2 and fz(x) I1 I2 corresponds to the following token sequence [+ COS mul + 2424 E-3 X1 X2 mul X1 x2]. Using this simple method, the sequence length scales linearly with the dimensionality of the system, i.e., the number of variables. While this is unproblematic for small dimensions such as D < Dmax = 6, it may impair the scalability of our approach." 4 As the encoder is only concerned with numeric input trajectories, its vocabulary only includes tokens for numbers.
Decoding strategy. At inference, we use beam sampling (Van Gael al., 2008) decode candidate equations, and select the candidate with highest reconstruction R² score. 5 For this, beam temperature an important parameter to control diversity- as the beam size increases, typically becomes useful to also increase the temperature to maintain diversity. Unless stated otherwise, we perform our experiments with a beam size of 50 and a temperature of 0.
Inference-time rescaling. During training, the model only observes trajectories integrated on the interval [1, 10] and starting from an initial condition sampled according to standard normal
A possible alternative would be to treat the decoding of each component function separate problem, adding specifier to the #08 (heginning of sequence) token to identify which component is be decoded.
BBeam search tends to produce candidates which all have the same skeleton. and only differ by small variations of the constants, leading to lack of diversity. Beam sampling ensures that randomness is added at each step of decoding leading to more diverse set of candidate expressions.
5

---

Published as a conference paper at ICLR 2024
distribution, described in Section 3. To accommodate for different scales of initial conditions and time ranges encountered during inference, we apply the affine transformations t + at + b to rescale the observed time range to [1, 10) and 2, + (t) 2.(0/2((ta) to rescale initial values unity. The prediction f that ODEFormer provides on inputs (f, x) are then transformed as .z.(,,] f: recover original units.
fi dt (2,(
dt
Optional parameter optimization. Most SR methods break the task into two subroutines (possibly alternating between the two): predicting the optimal "skeleton", i..., equation structure, and fitting of the constants contained in the skeleton. Just like Kamienny et al. (2022), our model is end-to-end in the sense that it handles both subroutines simultaneously. However, we also experimented with adding an extra parameter optimization step, as performed in methods such as ProGED (Omeje et al., 2023). We describe the details of the parameter optimization procedure in Appendix D) and denote this method as "ODEFormer (opt)".
EXPERIMENTS
Experimental tasks. For dynamical SR there is spectrum of reasonable tasks. We could simply be interested in finding some ODE f. whose solution approximates the observed trajectory zxt) on the observed time interval, even if f still differs from f on (parts of) their domain. A more ambitious goal closer to full identification is to find a f that also approximates the correct trajectories for unobserved initial conditions. This evaluation, which in our view is most meaningful to assess dynamical SR, is often absent in the literature, e.g., from Omejc et al. (2023). If an inferred f yields correct solutions for all initial values and time spans, we would consider as perfect identification of f. Accordingly, we evaluate the following aspects of performance in our experiments.
we compare the (noiseless, dense) ground truth trajectory with the trajectory obtained by integrating the predicted ODE from the same initial condition and on the same interval [1 as the ground truth.
Reconstruction:
Generalization: we integrate both the ground truth and the inferred ODE for new, different initial condition over the same interval [1, 7) and compare the obtained trajectories.
Metrics. When evaluating SR methods. the desired metric is whether the inferred ODE perfectly agrees symbolically with the ground truth expression f. However, such an evaluation is problematic primarily because (i) the ambiguity in representing mathematical expressions and non-determinism of "s iimpll i fy" operations modern computer algebra systems render comparisons on the symbolic level difficult, and (ii) expressions may differ symbolically, while evaluating to essentially the same outputs on all possible inputs (c.g., in the presence of a negligible extra term). Hence, comparing expressions numerically algebra systems include numerical evaluations in equality checks, which still require choosing a range on which to compare expressions (Meurer al., 2017). We consider the following metrics.
relevant range is more reliable and meaningful. Even modern computer
Accuracy: A classical performance metric for regression tasks is the coefficient of determination, (-x0, 1]. Since R2 is unbounded from below, average R across multiple predictions may be severely biased by even a single particularly poor outlier. We therefore instead report the percentage of predictions for which the R²2 exceeds threshold of 0.9 in the main text and show the distribution of scores in Appendix I.
defined as R² = E.(,, - -0)'
()) (v-
Complexity: we define the complexity of symbolic expression be the total number of operators, variables and constants. We acknowledge that this is a crude measure, for example, it would assign a lower complexity to exp(tan(z)) (complexit=33) than there is no agreed-upon meaning of the term "complexity" for a symbolic expression in this context so that simply counting "constituents" common in the literature, see e.g. the discussion in the benchmark by La Cava et al. (2021). Inference time: time to produce a prediction.
+ 2r (complexity=5). Yet.
Datasets. We evaluate ODEFormer on test set which is generated following the description in Section 3. Although the test set is sampled from the same distribution as the training data, the
presence of continuous constants within every sample (step 8 in the sampling procedure in Section 3) as well as i...d. initial values ensures that the test data does not overlap with the training data. In addition we also evaluate approach on two datasets not obtained via our generation procedure.

---

Published as a conference paper at ICLR 2024
Accuracy (R*x0.9)
80
60
40
20
2
0
2
6
20
40
60
50
100
150
200
dimensions
unary operators
# binary operators
points
reconstruction
clean data
noise, subsampling 5% noise, 25% subsampling
generalization
0% noise, 25% subsampling
Figure 3: Ablation study on synthetic data. We vary four parameters governing the difficulty of an example (from left to right): the dimension of the system, the total number of unary and binary operators, and the number of points the trajectory. In each panel, we average the results over dataset of 10,000 examples.
For this we first consider the "Strogatz dataset", included in the Penn Machine Learning Benchmark (PMLB) database (La Cava et al., 2021). consists of seven ODE systems and has been used as a benchmark by various SR methods, in particular those specialized
dynamical SR (Omejc et al., 2023). However, it has several limitations: (i) it is small (only seven unique ODEs, each integrated for 4 different initial conditions), (ii) it only contains 2-dimensional systems, (iii) it is not integrated with sufficient precision (Omeje et al., 2023), and (iv) its annotations are misleading (e.g., claiming that all systems develop chaos even though none of them does). Faced with these limitations. we introduce ODEBench, an extended benchmark curated from ODEs that have been used by Strogatz (2000) to model real-world phenomena as well as well-known systems from Wikipedia. We fix parameter values to obtain the behavior the models were developed for, and choose two initial conditions for each equation to evaluate generalization. ODEBench consists of 63 ODEs (ID: 23, 2D: 28, 3D: 10, 4D: 2), four of which exhibit chaotic behavior. We publicly release ODEBench with descriptions, sources of all equations, and well integrated solution trajectories - more details are in Appendix A.
Corruptions. For all datasets, also compare models on their robustness to two types of corruption: (i) to simulate imprecise observations we add measurements noise to the observations via zy (t,) - (1 + (f4) for j [D], E [N] and ~ (00 o); (ii) to simulate missing data and imperfect sampling processes, we drop a fraction p of the observations along the trajectory uniformly at random. We report results for various noise levels and subsampling fractions p. In addition, experiment when entire intervals of data are missing (e.g. due to sensor fault) can be found in Appendix F.
Baselines. In our experiments, extensively compare ODEFormer with (strong representatives of) existing methods described in Table 1. All baseline models are fitted to each equation separately, including separate hyperparameter optimization per equation. In contrast, we evaluate single, pre-trained ODEFormer model on all equations in all benchmarks. ODEFormer is hence the only model which does not require any hyperparameter tuning or prior knowledge on the set of operators to be used. Apart from ProGED and SINDy, all baselines were developed for functional SR. We use them for dynamical SR as described in Section 2, by computing temporal derivatives i; (t) via finite differences with hyperparameter search the approximation order and optional use of Savitzky-Sar Savgol filter for smoothing. For more details, please refer to Appendix E. Moreover, for performance comparison to the closely related univariate ODE inference model by Becker et al. (2023) see Appendix H.2.
5.1 RESULTS
Synthetic data We first assess how the performance of ODEFormer is affected by the dimensionality of the ODE system, the number of operators, and the number points the trajectory. Results on test set comprising 10,000 i.i.d. examples generated as described in Section are shown in Figure 3 we asses the effect the beam size in Appendix G. We make three observations:
Performance degrades with the system dimensionality as well as number of unary and binary operators as expected while ODEFormer is surprisingly insensitive to the number of points in the trajectory across the tested range.
Generalization accuracy is substantially lower than reconstruction accuracy expected, but at least for low owwdimensioonl systems we achieve non-trin con-trivial generalization (e.g., 60% generalization accuracy 85% reconstruction accuracy for ID).
ODEFormer copes well with subsampling, but suffers more from noisy trajectories. However, the effect generalization is smaller than that reconstruction.
7

---

Published as a conference paper at ICLR 2024
Table 1: Overview of models. f.d.: finite differences required, ode: method developed for dynamical
SR, T: munsformer-basede GP: genetic programming, MC: Monte Carlo, reg: regression
name
type
ode
f.d.
description
reference
ODEFormer
T
yes
n0
seq.to-seq. translation
ours
AFP
GP
no
yes
age-fitness Pareto optimization
(Schmidt & Lipson, 2011)
FE-AFP
GP
no
yes
AFP with co-evolved fitness estimates
(Schmidt & Lipson, 2011)
EHC
n1
yes
AFP with epigenetic hillclimbing
(la Cava, 2016)
EPLEX
GP
nO
yes
epsilon-lexicase selection
(La Cava et al., 2016b)
PySR
GP
nO
yes AutoML-Zero simulated annealing
(Crammer, 2023)
SINDY
reg
yes
yes
sparse linear regression
(Brunton et al., 2016)
FFX
reg
no
yes
pathwise regularized ElasticNet regression MC on probabilistic context free grammars
(McConaghy, 2011)
ProGED
MC
yes
yes
(Omeje et al., 2023)
Reconstruction results. We present results on both "Strogatz" and ODEBench in Figure 4. From top to bottom, we ranked methods by their average accuracy across all noise and subsampling levels. The ranking is similar on the two benchmarks and ODEFormer achieves the highest average score on both. The two leftmost panels show that ODEFormer is only occasionally outperformed by PySR when the data is very clean as noise and subsampling kick in, ODEFormer gains an increasingly large advantage over all other methods. In the two rightmost panels of Figure 4, we show the distributions of complexity and inference time. ODEFormer runs on the order of seconds, versus minutes for all other methods except SINDy, while maintaining relatively low and consistent equation complexity even at high noise levels. We show figures for all predictions in Appendix A.
aT 0.0
a-0,01
0-0.02
0-0.03
a 0.04
0-0.05
ODEFormer (opt)
p
p-
PySR
FFX
EHC
SINDY (poly)
SINDY (esc)
340
ProGED (poly)
ProGED
EPLEX
FE-AFP
SINDY
20
100
60
80 100
107
10
10a
107
10a
10º
10a
10
Accuracy (R² > 0.9)
% Accuracy (R2
0.9)
complexity
inference time [sec.]
(a) Reconstruction on Strogatz
ODEFormer (opt)
0.5
ODEFormer
PySR
FFX
SINDY (poly)
EHC
ProGED (poly)
+
ProGED
AFP
SINDy (esc)
FE-AFP
EPLEX
SINDY
Y!
20 40
60
80 100
20 40
60
80
100
10"
10?
10a
10a
10a
102
10a
% Accuracy (R² > 0.91
% Accuracy (R? 0.9)
complexity
inference time [sec.]
(b) Reconstruction on ODEBench
Figure 4: Performance comparison. We compare ODEFormer (and ODEFormer (opt) for additional parameter optimization) to existing methods, following the protocol described in Section 5. We present results for two values of the subsampling parameter p and six values of the noise parameter Whiskers box plot panels mark minimum and maximum values.
8

---

Published as a conference paper at ICLR 2024
Generalization results. We present general- ization results on ODEBench in Figure 5. Con- sistently across all models, accuracies drop by about half, meaning that half the correctly recon- structed ODEs do not match the ground truth symbolically. This highlights the importance of evaluating dynamical SR on different initial conditions. Note, however, that the overall rank- ings of the different methods rather consistent with the reconstruction results, and ODEFormer achieves the best results on average thanks to its robustness to noise and subsampling.
mx0.02
gH0.33
0=2.03
ODEFormer (opt)
p=0
10A
P 0.5
ODEFormer
PySR
ProGED (poly)
1/4
FFX
x
ProGED
146
AFP
FE-AFP
SINDy (esc)
EPLEX
SINDy (poly)
SINDY
HAY
25
50
25
50
% Accuracy (R² 0.9) % Accuracy IR² 0.9.
Figure 5: Generalization on ODEBench. We consider the same setting as Figure 4.
6 CONCLUSION
We foresee real-world applications of ODEFormer across the sciences, for hypothesis generation of dynamical laws underlying experimental observations. However, in the following we also highlight several limitations of the current method, opening up interesting directions for future work.
First, we only considered first order ODEs. While any higher-order ODE can be written as a system of first order ODEs, this does not immediately allow ODEFormer to make predictions based only on a solution trajectory, since we would still need to approximate time derivatives up to the order of the ODE. While possible principle via finite differencing schemes, we would suffer similar drawbacks as other methods that rely on finite differences from noisy, irregularly sampled data.
Second, ODEFormer only works when all variables are observed. In real-world settings, some relevant variables may be unknown or unobservable. For example, inferring chemical kinematics may be challenged by the difficulty in measuring the concentration of reaction intermediates (Burés & Larrosa, 2023). We may circumvent this issue by randomly masking variables during training (replacing their observations by a dedicated token), to emulate unobserved variables. While this raises questions around identifiability and robustness, we plan to explore masking in future work. This technique could also handle higher-order ODEs by considering the derivatives as unobserved.
Third, based the four representatives in ODEBench, ODEFormer (as well as all other benchmarked models) struggles with chaotic systems, which have been argued to provide good assessments for data- driven modeling and forecasting (Gilpin, 2021). For chaotic systems, understanding properties of the attractor is often more desirable: nearby trajectories diverge exponentially, rendering identification notoriously challenging. Dynamical law learning from short (transient) trajectories for chaotic systems remains an interesting direction for future work.
Lastly, all existing methods for dynamical SR, including ours, perform inference based on single observed trajectory. In our opinion, one of the most promising directions for future work is to enable inference from multiple solution trajectories of the same ODE. Key benefits may include averaging out possible noise sources as well as improving identifiability, as we "explore the domain of f". However, initial experiments with various forms of logit aggregation in ODEFormer's decoder during inference did not yield convincing results. In future work, we plan to exploit a cross-attention encoder to combine the embeddings of the different trajectories to leverage the combined information, the spirit of Liu et al. (2023).
In addition to developing solutions to these limitation another exciting direction for future work concerns the extension of the framework towards other classes of differential equations, e.g. partial or stochastic DEs. Although conceptually straightforwards, finding scalable data generation procedures for these classes may be challenging yet would considerably expand the set of possible applications.
We conclude by emphasizing that the difficulty and potential ambiguity in dynamical law learning calls for caution during deployment. Methods like ODËFormr primarily serve as hypothesis generators, ultimately requiring further experimental verification, something that can be done well based on ODEs (rather than black-box models).

---

Published as a conference paper at ICLR 2024
7 REPRODUCIBILITY STATEMENT
The reproducibility of our work ensured through several means. First and foremost, all code, model weights, and created benchmark datasets will be made publicly available http://ituub.comsssassolii/////e
together with notebooks to directly reproduce key results, as well pip-installable package for easy usage. We also describe detail the data generation in Section 3, our architecture and model choices in Section 4, and additional details about training in Appendices and D.
ACKNOWLEDGEMENT
We thank Luca Biggio, François Charton and Tommaso Bendinelli for insightful conversations. We acknowledge funding from the EPFL's AI4science program (SA), Helmholtz Association under the joint research school "Munich School for Data Science MUDS" (SB), Swiss National Science Foundation grant (310030 21002222222 212516) (AM), NCCR Catalysis (grant number 180544), and a National Centre of Competence in Research funded by the Swiss National Science Foundation (PS). The authors gratefully acknowledge the Gauss Centre for Supercomputing e. V. (www.gauss-cer gauss-centre.eu) for funding this project (application number 29063, project "dynadis") by providing computing time on the GCS Supercomputer JUWELS (Jülich Supercomputing Centre, 2021) at Jülich Supercomputing
Centre (JSC).
REFERENCES
Nikos Aréchiga, Francine Chen, Yan-Ying Chen, Yanxia Zhang, Rumen Iliev, Heishiro Toyoda, and Kent Lyons. Accelerating understanding of scientific experiments with end to end symbolic regression. ArXiv, abs/2112.04023 2021.
Karl Johan Äström and Peter Eykhoff. System identification-a survey. Aulomatica, 722):123-162 1971.
Steven Atkinson, Waad Subber, Liping Wang, Genghis Khan, Philippe Hawi, and Roger Ghanem. Data-driven discovery of free-form governing differential equations. arXiv preprint 1rX1v:1910.05117, 2019.
Sören Becker, Michal Klein, Alexander Neitz, Giambattista Parascandolo, and Niki Kilbertus. Predicting ordinary differential equations with transformers. In International Conference Machine Learning, pp. 1978-2002. PMLR, 2023.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In International Conference on Machine Learning, pp. 936-945. PMLR, 2021.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 1 113(15):3332-39)
2016. ISSN 0027-8424. doi: 1..1073/pna.15117388 URL
https: www. phas..n........... org/ wwntent//13/15///1
Jordi Burés and Igor Larrosa. Organic reaction mechanism classification using machine learning. Nature, 613(7945):689-665 2023.
Anja Butter, Tilman Plehn, Nathalie Soybelman, and Johann Brehmer. Back to the edition. arXiv preprint arXiv:2109.1044 2021.
Francois Charton. Linear algebra with transformers. Transactions on Machine Learning Research, 2022.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 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.
10

---

Published as a conference paper at ICLR 2024
Miles Cranmer. Interpretable machine learning for science with pysr and symbolicregression, jl. arXiv preprint arXiv:2305.01535 2023.
Stéphane d' Ascoli, Fierre-Alexandi 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 Proceedings of Machine Learning Research, pp. PMIR, 17-23 Jul 2022. URL https: /twpe/eoooeeeeeeeoeooooooooooooc ings ngs. wwm.m.rresss... v1 62/ d- d-ascoli224, html.
A Davies, P Velickovic, L. Buesing, S Blackwell, D Zheng. N Tomasev, Tanburn, P Battaglia, C Blundell, Juhasz, et al. Advancing mathematics by guiding human intuition with ai. Nature, 2021.
Brian de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, J Nathan Kutz, and Steven Brunton. PySINDy: Python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49): 2020.
Jonas Degrave, Federico Felici, Jonas Buchli, Michael Neunert, Brendan Tracey, Francesco Carpanese, Timo Ewalds, Roland Hafner, Abbas Abdolmaleki, Diego de Las Casas, et al. Magnetic control of tokamak plasmas through deep reinforcement learning. Nature, 602(7897):414 419, 2022.
John R Dormand and Peter J Prince. A family embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26, 1980.
X Duan, JE Rubin, and D Swigon. Identification of affine dynamical systems from single trajectory. Inverse Problems, 36(8):085004, 2020.
Rudresh Dwivedi, Devam Dave, Het Naik, Smiti Singhal, Rana Omer, Pankesh Patel, Bin Qian, Zhenyu Wen, Tejal Shah, Graham Morgan, et al. Explainable ai (xai): Core ideas, techniques, and solutions. ACM Computing Surveys, 55(9):1-33, 2023.
Stefan Elfwing, Eiji Uchibe, and Kenji Doya. Sggmoid-weighted linear units for neural network function approximation in reinforcement learning. Neural networks, 107:3-11, 2018.
Sébastien Gaucel, Maarten Keijzer, Evelyne Lutton, and Alberto Tonda. Learning dynamical systems using standard symbolic regression. In Genetic Programming: 17th European Conference, EuroGP 2014, Granada, Spain, April 23-25, 2014, Revised Selected Papers 17, pp. 25-36. Springer, 2014.
William Gilpin. Chaos as interpretable benchmark for forecasting and data-driven modelling. In fiitty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2), 2021.
Roger Guimerà, Ignasi Reichardt, Antoni Aguilar-Mogas, Francesco Massucci, Manuel Miranda, Jordi Pallarès, and Marta Sales-Pardo. A bayesian machine scientist to aid in the solution of challenging scientific problems. Science advances, 6(5):caav6971, 2020.
Christopher Hahn, Frederik Schmitt, Jens Kreber, Markus N Rabe, and Bernd Finkbeiner. Teaching temporal logies neural networks. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021. URL https:// openreview net/ forum 1d=dd=okkff4505 f4byz.
Frédéric Hamelin, Abderrahman Iggidr, Alain Rapaport, and Gauthier Sallet. Observability, identifia- bility and epidemiology a survey, arXiv preprint arXiv:2011.12202 2020.
A.C. Hindmarsh and Lawrence Livermore Laboratory. ODEPACK, a Systematized Collection of ODE Solvers. Lawrence Livermore National Laboratory, 1982. URL https:: h/op:ss///okkss// google. de/books ?id= 9XWPHWEACAAJ.
Jülich Supercomputing Centre. JUWELS Cluster and Booster: Exascale Pathfinder with Modular Supercomputing Architecture at Juelich Supercomputing Centre. Journal of large-scale research fa- cilities. 7(A188), 2021. doi: 10.17815/jlsff-7 7813. URL http:/ dxx doi. .or. 110g//0...00 17815/ 71
11

---

Published as a conference paper at ICLR 2024
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn /unyaauvvnakkol,, Russ Bates, Augustin Zidek, Anna Potapenko, et al. Highly accurate protein structure prediction with alphafold. Nature, 596(78733:583-558) 2021.
Fierre-Alexande Kamienny, Stéphane d' Asooll, Guillaume Lample, and Francois Charton. End-to- end symbolic regression with transformers. In Alice Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022. URL https: /pwpneeeeeeeeeeeeeeeeeeeeeeeeeee net/ / forum? God.dddoGGoHHGH......
George Em Karniadakis, loannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 2015.
Michael Kommenda, Bogdan Burlaeu, Gabriel Kronberger, and Michael Affenzeller. Parameter identification for symbolic regression using nonlinear least squares. Genetic P Programming and Evolvable Machines, 21(3):41--50, 2020. doi: 10. 10.10007/10711-0010-00777
https:
1000 010-09091199000000 09371-3.
Gabriel Kronherger, Lukas Kammerer, and Michael Kommenda. Identification of dynamical systems using symbolic regression. In Roberto Moreno-Díaa Díaz, Franz Pichler, and Alexis Quesada-Arencinn (eds.), Computer Aided Systems Theory - EUROCAST 2019, pp. 370-377, Cham, 2020. Springer International Publishing. ISBN 978-3-030-450933
William La Cava, Kourosh Danai, and Lee Spector. Inference of compact nonlinear dynamic models by epigenetic local search. Engineering Applications of Artificial Intelligence, 55:292-306, 2016a.
William La Cava, Lee Spector, and Kourosh Danai. Epsilon-lexicase selection for regression. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 741-748, 2016b.
William La Cava, Tilak Raj Singh, James Taggart, Srinivas Suri, and Jason H Moore. Learning concise representations for regression by evolving networks of trees. In International Conference on Learning Representations, 2018.
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabricio Olivetti de Franca, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason Moore. Contemporary symbolic regression methods and their relative performance. Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmurks Track (Round 1), 2021.
William G La Cava. Automatic Development and Adaptation of Concise Nonlinear Models for System Identification. PhD thesis, University of Massachusetts Amherst, 2016.
Guillaume Lample and François Charton. Deep learning for symbolic mathematics. In International Conference on Learning Representations, 2019.
Mikel Landajuela, Chak Shing Lee, Jiachen Yang, Ruben Glatt, Claudio P Santiago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. unified framework for deep symbolic regression. Advances in Neural Information Processing Systems, 1..33985-33998 2022.
Yuchen Liu, Natasha Ong, Kaiyan Peng, Bo Xiong, Qifan Wang, Rui Hou, Madian Khabsa, Kaiyue Yang, David Liu, Donald S Williamson, et al. Mmvit: Multiscale multiview vision transformers. arXiv preprint 000000230,.00104 2023.
Nour Makke and Sanjay Chawla. Interpretable scientific discovery with symbolic regression: review. arXiv preprint arXiv:2211.10873, 2022.
Trent McConaghy. Ffx: Fast, scalable, deterministic symbolic regression technology. In Genetic Programming Theory and Practice IX, pp. 230-260. Springer, 2011.
Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondrej Čertík, Sergey B Kirpichev. Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. Sympy: symbolic computing in python. PeerJ Computer Science, 3:3103, 2017.
12

---

Published as a conference paper at ICLR 2024
Hongyu Miao, Xiaohua Xia, Alan S Perelson, and Hulin Wu. On identifiability of nonlinear ode models and applications in viral dynamics. SIAM review, 53(1):3-39, 2011.
Terrell Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio Santiago, Daniel Faissol, and Brenden K Petersen. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. Advances in Neural Information Processing Systems, 34:24912- 14224922-24923, 2021.
Jorge Nocedal and Stephen Wright. Numerical Optimization. Springer, New York, NY, USA, 2e edition, 2006.
Nina Omeje, Boštjan Gec, Jure Brence, Ljupčo Todorovski, and Sašo Džeroski. Probabilistic grammars for modeling dynamical systems from coarse, noisy, and partial data. Researchsquare preprint r5-2678362, 2023. doi: https://doi.org/ 3ttp://0..11000/3....22222535550
Fabian Pedregosa, Gaël Varoquaux. Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jak Vanderplas, Alexan- dre Passos, David Cournapeau, Matthicu Brucher, Matthicu Perrot, and Édouard Duchesnay. Scikit-learn Machine learning in Python. JMLR, 2011.
Brenden K Petersen, Mikel Landajuela Larma, Terrell N Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathemat- ical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2020.
Lena Podina, Brydon Eastman, and Mohammad Kohandel. Universal physics-informed neural networks: Symbolic differential operator discovery with sparse data. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Pro- ceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 27948-27955 PMLR, 23-29 Jul 2023.
Stanislas Polu and llya Sutskever. Generative language modeling for automated theorem proving. arXiv preprint 0tt::22009033933 2020.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-code: Discovering closed-foor foom odes from observed trajectories. In International Conference Learning Representations, 2021.
Xing Qiu, Tao Xu, Babak Soltanalizadeh, and Hulin Wu. Identifiability analysis of linear ordinary differential equation systems with a single trajectory. Applied Mathematics and Computation, 430: 127260, 2022.
Markus Quade, Markus Abel, Kamran Shafi, Robert K Niven, and Bernd R Noack. Prediction of dynamical systems by symbolic regression. Physical Review E, 94(1):012214, 2016.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informer neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Michael Schmidt and Hod Lipson. fiteffitess pareto optimization. In Genetic programming theory and practice VIII. pp. 129-146. Springer, 2011.
Philipp Scholl, Aras Bacho, Holger Boche, and Gitta Kutyniok. The uniqueness problem of physical law learning. In ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1-5. IEEE, 2023.
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, 58877839):604-669, 2020.
Shelby Stanhope, Jonathan E Rubin, and David Swigon. Identifiability of linear and linear-in- parameters dynamical systems from single trajectory. SIAM Journal on Applied Dynamical Systems, 13(4): 172--1815, 2014.
Steven H. Strogatz. Nonlinear Dynamics and Chaos: With to Physics, Biology, Chemistry and Engineering. Westview Press, 2000.
13

---

Published as a conference paper at ICLR 2024
Erico Tjoa and Cuntai Guan. A survey explainable artificial intelligence (xai): Toward medical xai. IEEE transactions on neural networks and learning systems, 32(11):4793-4813 2020.
Tony Tohme, Dehong Liu, and Kamal Youcef-Toumi. Gsr: A generalized symbolic regression approach. Transactions Machine Learning Research, 2022.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: physics-inspired method for symbolic regression. Science Advances, 6(16):eaay2631, 2020.
Silviu-Marian Udrescu and Max Tegmark. Symbolic pregression: Discovering physical laws from raw distorted vidco. Physical review. E, 103 4-1:043307, 2021.
Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. Symboliegpt: A generative trans- former model for symbolic regression. arXiv preprini arXiv:2106. 14131, 2021.
Jurgen Van Gael, Yunus Saatci, Yee Whye Teh, and Zoubin Ghahramani. Beam sampling for the infinite hidden markov model. In Proceedings of the 25th international conference on Machine learning, pp. 1088-1095, 2008
Martin Vastl, Jonáš Kulhánek, Jirí Kubalík, Erik Derner, and Robert Babuška. Symformer: End-to- end symbolic regression using seaeformer-bbsee architecture. arXiv preprint arXiv: 2022.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in neural information processing systems, pp. 5998-6008, 2017.
Alejandro F Villaverde, Antonio Barreiro, and Antonis Papachristodoullu. Structural identifiability of dynamic systems biology models. PLoS computational biology, 12(10))e1005153, 2016.
Marco Virgolin, Tanja Alderliesten, Cees Witteveen, and Peter AN Bosman. Improving model-based genetic programming for symbolic regression of small expressions. Evolutionary computation, 29 (2):211-22) 1-237, 2021.
Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, al. Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods, 17(3):66--22,, 2020.
Juliane Weilbach, Sebastian Gerwinn, Christian Weilbach, and Melih Kandemir. Inferring the structure of ordinary differential equations. arXiv preprint arXiv:2107.07345 2021.
14

---

Published as a conference paper at ICLR 2024
A ODEBENCH
ODEBench features a selection of ordinary differential equations primarily from Steven Strogatz's book "Nonlinear Dynamics and Chaos" with manually chosen parameter values and initial condi- tions (Strogatz, 2000). Some other famous known systems have been selected from other sources such as Wikipedia, which are included the dictionary entries well. We selected ODEs primarily based on whether they have actually been suggested as models for real-world phenomena as well as on whether they are 'iconic' ODEs in the sense that they are often used as examples in textbooks and/or have recognizable names. Whenever there were 'realistic' parameter values suggested, we chose those.
In this benchmark, we typically include only one set of parameter values per equation. Many of the ODEs in Strogatz' book are analyzed terms of the different limit behavior for different parameter settings. For some systems that exhibit wildly different behavior for different parameter settings, include multiple sets parameter values as separate equations (e.g.. the Lorenz system in chaotic and non-chaotic regimes). For each equation, we include two sets of manually chosen initial conditions.
There are 23 equations in one dimension, 28 equations two dimensions, 10 equations three dimensions, and 2 equations in four dimensions. This results in total of 63 equations, 4 of which display chaotic behavior. We provide the analytical expressions and initial conditions in Tables to 4, visualizations of single trajectory for each ODE in Figure 6, and ODEFormer's predictions for each ODE in Figure?. 7.
15

---

Published as a conference paper at ICLR 2024
ID
System description
Equation
Pammeters
Tnitaa values
RC-sirveit (churging cupaciion)
0.7. 231
[100]. [3.54]
Population growth (nuive)
0.23
[4.78], (0.87]
Population growth with carrying pacily
(1
0.79. 74.3
[7.3]. (2100)
BC circuit with non linear resistoe (charging eapecitor)
-0.5 +
20
0.96
[0.8],
60 " +1
Velocity of falling object with air resistance
9.81, 0.0021051115
[7.0)
Autocatalyıi with one fixed ahun- dunt chemical
9120
2.1,0.5
[0.23]
Compertz low for tomor growih
coing log (cymg)
0.032
[1.33]
Logistic equation with Allee effect
coitg -1+ 50
0.14, 130.0,
(2.1)
Language death model for two lan-
Go c1 20) 6139
0.32 0.28
[0.4]] (0.55]
guages
language death
romgg (1 20)
2g 20-00111 (0) =2
0.2,
[0.3]]
two languges
slowing
-xg
[34]
tical mechanies)
(simple)
Godgg 0120
[1..)]
Overdamped bead on a rotating
6o (e) cos (#0) sijin (#o)
0.081, 9.7
[33]]
hoop
Rudworm outhroak model with pro- dation
rging (1 %
5
0.78, 81.0, 21.2, 21.0 09
[2.76]. (23.3)
15
Budworm outbresk with predution (dimendionees))
collig
(1 -
0.4, 95.0
[44.3] [4.5]
16
Landuu equation (ypical lime scile =
cortig 01
0.1. 0.000
[0.94] [1.65]
Logistic equation with harvest- ingitsshing
(1
0.4, 100.0, 0.3
[14.3]
Improovd logistic equation with vestingllishing
coltig
0.4, 0.34,
[21 [44.1]
Improved logistic equation with vestingttshingg (dimensionless)
2010
20
20)
0.08, 0.8
[0.33] [0.03]
x 60 01
Autocatalytic gene switching (di-
0.1, 0.55
[0002], [0.25]
mensionless)
#0+1
Dimensionallyy reduced SIR infec- tion dead people (dimen- sionless)
60 012g
#0
[00],
protein ex- pression (positive leodback, hasal promooter expression)
Hysteretic
6o
123.0,
16.3]
Overdamped with loepe/onneliees/oephoo tion (dimensinnless)
pendalum driving
60 (stg)
constans
june-
Table 2: Scalar ODEs in ODEBench.
16

---

Published as a conference paper at ICLR 2024
ID
System description
Equation
Purmeeters
Enitall values
24
Harmonic oscillator without damp-
21
104, 4002), 10.0, 0.1
ing
25
Harmonic oscillator with dhmping
#1
4.4,0.43
[0.12,00443) (0.0, -0.3
-00itg - c1.r
26
Volterr competition model sheeps and
#0 (eg 0141 - =20)
30, 2.0, 2.0
[50, (2.3,,)))
(Sirogatz
#1 (#g 20 23)
rabbits)
2p (eg 21)))
27
Lotks- Volterr simple Wikipedia)
tos
on
c3
1.54. 1.45, 3.0, 1.62
[8.3. 3.4,, [0.4, D.65)
28
Pendulum without friction
[-19, 0.0), [0.3.
-cn sin (xg)
29
Dipole fixed point
ro
=ag + of
0.65
[3.2 (1330
#0 (-coaggt + 21)
30
RNA molecules catalyzinng each oth- replication
1.61
0.04). (0.1, 0.21]
21 (-Co202) 2o)
model only honlthy
-00302]
[7.2, 0.9)) (200, 12.4]
and sick
co2g21
32
Dumped double well escillator
F1
0.18
[-L8, (58.
-m0m
40
+
-00
(z1)
33
Cliider (dimensionless
0.08
[50, 0.0) (9.81,
20
2g
34
Friztionless bead rotating hoop
0.33
(-1.2, -0.1
(limencionless)
(-r0 + coa (#o)) (mp)
(20) (#r)
35
Rotational dynamics of an object in shear flow
(co sin" (x1) eos? (x1) (#r)
4.2
[1113, 0.2] (2:4, 0.21 1.7]
damping.
[0,45, (13) 0.8]
no driving (dimensivnless)
- ooe (spp) 21 sin (mo)
21
37
Van der Pol oscillatoe (standard form)
(=) 11
0.43
(.0 (0.1,
38
der oscillator (impliifee form from Stragat))
3
3.37
007,000 (-1.1 -07]
39
Glycolytic oscillator, e.s.. ADP and IeP in yeast (dimensionless)
cpz1
x8x= 20
2.00
(..1 0.71
40
Durfing qquation (weakly non-lineze oscilliation)
21
6021
(1- =1)
(063, [0.2,03)
41 Cell cycle model by Tyson for teraction protein edo2 and cyclin (dimensionless)
to
(a1 + =0) (-) (-n0 11) -2g
15.1 0,001,.0
0.38 10.02, 1.2)
62 - 20
42 Reduced model for chlorine diotide iodine malenic rection sionless)
89,
0.35).35) 7.8)
45 Driven pandulum with lincar damp- ing / Josephson junction climension-
21
1,67, 0.64
[1.47, [1.47.02]] (-10).005))
6o
c1
(x0)
pesdulio with quadratic dumping (dimensioness))
1147,-02j.1-19. 0.03)
c1 [xi] sin
60- zo) mo27
45
Isothermal autocatalytic reaction medel Gray and Scot 1985 (di- metsionless)
0.5.0.02
[.A. (0.32, (0.32.0.44 0.64]
co (zo #1) (rO)
46
Interacting bur magnets
0.33
[0.56 [0:43, 1,21]
sin (ig 21 sin (#11)
47
Hinneular divley model oscilla- tions)
c0 11 -21
4.8%
Jo65, 10.3)
20
48
Dacterial respiration model for nu- trients oxygen levels
Co
c147 11
Ro
0.., 11.25
[0.1, 11332,2211
02
6187+1
0187 20 (eg +144
49
Brussclatoe: hypothetical chemical oscillatioon medel dimessionless)
c1
3,0,3.1 3.1
-07,114)
Chemical oscillator model
+
[0.14, 1..55..66
Schnackenberg 1999 (dimassion-
6 *6"1
less)
co + sin (#1 coe (x))
51 Oscillator death model by Emnen Kapell 1990
1432,0.972
(22)0066) 0.12)
e1 + (111) ooe (ing)
Table 3: dimensional ODEs in ODEBench.
17

---

Published as a conference paper at ICLR 2024
description
Equation
Initial
<a (-20 + x1)
Maxwell-Bie Bloch equations tlaser dy- namics)
(2gag
0.1,
11.3. 0.89), 1.1
(-c32021
+
c420
eg 1g
Model for apoptosis (cell death)
0122 (oa #1)
7,95,
10,005, 0.26
+22)
1andr
0.0027]
Lurenz well behaved periodic regime
9120 3002
5.1, 120. 1.67
[2.3. 8.I. 12.4). (10.0.
=c222 #.202
2000, 30.0
eg (=#0 +-21)
Lorenz equations complex, peri- udio regimme
41 2022
12.3. 12.41,
20.0, 30.00
(-20 #1)
Lorenz equations standard parame-
91 2022 -
10.0, 28.0, 8/3
12.3. 8.1. 12.4), (10.0. 20.0, 30.0)
ters (chaotac)
-0289 + 2g
(-21
57
Rüssler atractor (stable fixed point)
(cg21 + 26)
0.2. 0.2, 5.3.50
12.3. 11, [.1. 4.0.1. 2.11
cy (c) #.g (-0g (=0)+ mg))
<g [-21 22)
Rössler (periodic)
cg +0)
0.2, 5.7,
12.3. O.L.0.8.0.....00
(-02 [0] +.2 (-02 xg))
2.1]
<g [-#1 22)
Rössler (chantic)
4g [0g2) + x6)
5.7.
[2.3 4-0.1,411
[ci
sq))
2.1]
-caa1 + (-0 zg]
Aizawn attracice (chaotic)
<g 21 (-e1 + #2)
0.95. 0.7 0.65, 3.5,
[0.1. 0,05, 0.55, [-0.3,
6g22 0 egageg 1/3a - (+0 +2)
(ca22 +
61 Chen Lee uitractor; system for gyu motion with feclback control of body
4121 + 2032
5.0, 10.0, -38, 3.0
[15. 15, 155, 18. 14,
6222 1971
62 Dinocular risulry model with slup- tation (oscillations)
<3 (it m1]
0.89, 0.9.0.4. L44 1.0
12.25, -0,5, -1.13, 044).
-22 210
10.342, -0.451,
63 (12 F3)
-c; Mg22
-op01 + 218p22 SEIR infection model (proportios)
0.47, 0.28, 0.3
10.6, 0.3, 0.09 0.1). [0.4, 0.3, 0.25,
0222
Table 4: and dimensional ODEs in ODEBench.
18

---

Published as a conference paper at ICLR 2024
Orgar
Loplte
www.u
Butwer www.eeewettotcccoooooooooo pediiie
Landa tppeiitelotttee
Lapakogptto vith
Iepanddingndoo mih
Tumanay 3694
ARt00 preadkiing
Elottined bedked peglo
X--
Crmiangi reddad star
Car cralin
Langlng condedt
Cale Iedride and
and and sotue
de
Mmmm
Figure 6: Solution trajectories of all equations in ODEBench for one of the initial conditions.
19

---

Published as a conference paper at ICLR 2024
Figure 7: Predictions of ODEFormer for all equations in ODEBench for the first set of initial conditions.
20

---

Published as a conference paper at ICLR 2024
B IDENTIFIABILITY
Traditionally, when inferring dynamical laws as ODEs from observations, one assumed a parametric form of the ODE to be known. In addition, typical assumptions included that the state of the system may not be fully observed, and that one may not be interested in identifying the system fully (all parameters), but that there are only certain (combinations of) parameters that need to be identified. For example, in typical epidemiological model of disease spread, e.g., modelling the fractions of susceptible, infected, and recovered people (SIR model) in population, one may only be able to measure the number of infected people, but may also only be interested in estimating the reproduction rate (a ratio of two parameters in the full equations). Hence a typical identifiability question would start from known system, known observation function, and a specific target quantity. Researchers have then developed various general methods and procedures to decide whether such a query is solvable, i.e., whether the target is identified from the observations within the assumptions (Aström Eykhoff, 1971; Miao et al., 201 I; Villaverde et al., 2016; Hamelin et al., 2020).
The study of identifiability from observations within larger non- non-parametrie function classes from full state observations has only been taken up more recently. For example, linear autonomous systems as well as autonomous systems that are linear in parameters are fairly well understood (Stanhope et al., 2014; Duan et al., 2020) with Qiu et al. (2022) recently essentially providing closure to this question. A broad summary of these findings is that linear (in parameters) autonomous systems are almost always (almost surely) identifiable from single trajectory for many reasonable measures (probability distributions) over parameters, i.e., over ODE systems.
When it comes to non-parametric classes such analytie, algebraic, continuous, or smooth functions, Scholl et al. (2023) have recently presented the first detailed analysis and results for a broad class of scalar PDEs. These include ODEs special cases, however their results only apply to scalar ODEs. Even though they consider ODEs beyond first order, this still does not include multivariate ODE systems. In the scalar case, identification from a single solution trajectory is possible for analytic functions f, but essentially impossible for continuous (or smooth) f (Scholl et al., 2023).
The function class considered by ODEFormer, i.e., the distribution implied by our generator, pre- dominantly contains real analytic functions. but not exclusively since + 1/2 is analytic only on \AO {0). Hence, it is not clear in which category we fall regarding identifiability.
Crucially, all these theoretical results assume the entire continuous, non-noisy solution trajectory to be known. Little known for discrete and noisy observations, where identifiability likely turns from yes/no question into one of probabilistic claims given prior over functions and the assumed noise model. Hence, current theoretical results do not conclude whether non-linear ODEs can practically be inferred from data.
C TRAINING DETAILS FOR ODEFORMER
We optimize the cross-entropy loss between the predicted sequence of mathematical symbols and the target sequence. During training, (the night-hand side of) ODEs are represented in prefix format instead of the more familiar infix format in order to relieve the model from having to predict correct opening and closing parenthesis (as an example: "3 (xo + x1)" corresponds to ", +, Tg 31 in prefix notation). We note that different target sequences may be mathematically equivalent but do not attempt to resolve this ambiguity on the training data set as it does not hamper model optimization in practice.
We use the Adam optimizer (with default parameters suggested by Kingma & Ba (2015)), with a learning rate warming up from 10-7 to 2 10-4 across the initial 10,000 steps and subsequent decaying governed by a cosine schedule for the next 300,000 steps. The annealing cycle then restarts with damping factor of 3/2 as per Kingma & Ba (2015), resulting in approximately 800,000 optimization steps. We do not use any regularization such as weight decay or dropout. To efficiently manage the greatly varying input sequence lengths, we group examples of similar lengths in batches, with the constraint that each batch contains 10,000 tokens. Our model is trained on set of about
'The intuition here is that linear, polynomial, or more broadly even real analytic scalar univariate functions can be uniquely extrapolated to all of RR when they are known on any open infinitely many ways extrapolate continuous smooth functions heyond any interval.
interval. On the other hand, there are
21

---

Published as a conference paper at ICLR 2024
50M examples pre-generated with 80 CPU cores. When single NVIDIA A GPU with 80GB memory and 8 CPU cores, ODEFormer's 's training process takes roughly three days.
To enforce model robustness to signal corruptions contaminate training trajectories x(ti) with noise according to (t,( + (1+ + (tz)(t) for E [D].i E [N] and Ç ~ NN
An overview of the model architecture can be found in Table 6.
D OPTIONAL CONSTANT OPTIMIZATION IN ODEFORMER (OPT)
In contrast to all baseline models, ODEFormer is a pretrained model and predicted ODEs are not explicitly fit to the data observed at inference time. However, similar to Kamienny et al. (2022) we can post-hoc optimize the parameters of a predicted ODE to improve the data fit. Although parameter estimation for dynamical system is known be challenging inference problem, algorithm (BFGS) (Nocedal & Wright, 2006) as implemented in s scipy- . imi ze mi nimi (Virtanen et al., 2020) and thus opt for comparatively simple local, gradient-basee method in the hope that the parameter values predicted by ODEFormer only need slight refinement. The optimizer solves the following problem
the Fooyedn--
arg min loss( (30(4)). -.. (1 solve solve_ivp( (2; 0), zo s<lo), ('o ; ,1,))
where {pI,.... pR) denotes the set of parameters of the ODE f that was predicted by ODEFormer, and where (z(44) ..... .x())) represents the (potentially noisy) observations. We use the negative werihtce-weighted R²² score as optimization loss.
.(.))(()))))
E EVALUATION OF BASELINE MODELS
Hyperparameter optimization. All baseline models are fitted to each trajectory separately and each fit involves a separate hyperparameter optimization. Hyperparameters that are searched over listed in Table 5, all other hyperparameters are set their respective default values. For each combination of hyperparameters, the model is fitted on the first 70% and scored on the remaining 30% of trajectory. To reduce runtime, we parallelize optimization according to Gr dSea rchCV from scikit- Malkit-learz learn (Pedregosa et al., 201 I) and set the number of parallel jobs to min (# combinations, # cpu cores 48)). After selecting the combination with highest Rt score, the final model is fitted on the full trajectory.
Finite difference approximations. Except for ProGED, all baseline models require approxima- tions of temporal derivatives of all state variables of an ODE system as regression targets. To estimate temporal derivatives we use the central finite difference algorithm as implemented by FiniteDi fference in the pysindy software package (de Silva et al., 2020) and include the approximation order in the hyperparameter search. For a fair comparison on noisy trajectories we extend the hyperparameter search to also include optional smoothing of trajectories with a Savitzky- Savgol filter with a window length of 15 as implemented by Smoot iiweeedddiiieiiiii ference (de Silva et al., 2020).
Vector-valued functions. Some of the baseline implementations (AFP, FE-AFP, EPLEX, EHC, FFX) do not readily support vector-valued functions ( RDD D) but only scalar-valued functions (f D + IR.. To evaluate these baselines on systems of ODEs, we run them separately for each component fi R + R2 of the system and combine the predictions for all components i E (I...., D} via the Cartesian product (st,.. ...., _,,k] [ (S....., [KK where K, represents the number of predictions, e.g., the length the Pareto front, obtained for component i.
Candidate selection. In symbolic regression, one typically faces a trade-off between accuracy (how well the function//rajectorry is recovered) and complexity of the proposed expression. There are different strategies in the literature to select a single, final equation from the accuracy-complextty Pareto front, which may bias comparisons across methods along one or the other dimension. For fair comparison,
evaluate all equations of model's Pareto front and pick the final equation based on R² score between (potentially noisy) observed trajectory and the integrated solution of the predicted ODE.
22

---

Published as a conference paper at ICLR 2024
Table 5: Hyperparameter names and values for optimization of baseline models. For FFX and PySR, we optimize over finite difference order and smoother window length but no additional hyper parameters.
Model
Hyperparameter
Values
finite difference order smoother window length population size generations
2, 3, 4
All models
None, 15
100, 500, 1000
AFP
2500, 500, 250
[n, V. +, *, /, exp, log. 2, 3, sqrt],
operators
[n. V. &, /, exp. log, 2, 3, sqrt, sin, cos]
population size generations
100, 500, 1000
EHC
1000, 200, 100
[n, V, +, *, /. exp, log. 2, 3, sqrt].
operators
[n, v. 4. * exp, log, 2, 3, sqrt, sin, cos]
population size generations
100, 500, 1000
EPLEX
2500, 500, 250
[n. -, *, 1, exp, log, 2, 3, sqrt],
operators
[n, v. +, $, 8,/, exp, log, 2, 3, sqrt, sin, cos]
population size generations
100, 500, 1000
FE-AFP
2500, 500, 250
[n. v. * exp. log, 3. sqrt],
operators
v. ..4. exp. log. 2. 3, sqrt, sin, cos]
universal, rational simplerational, trigonometrie polynomial
ProGED
grammar
polynomial degree
2, 3, 4, 7, 8, 9, 10 [polynomials],
SINDy
basis functions
[polynomials, sin, cos, expl.
[polynomials, sin, cos, exp, log, sqrt, 1/x]
optimizer threshold
0.05, 0. 1, 0.15
optimizer alpha
0.025, 0.05, 0.075
optimizer max iterations
20, 100
polynomial degree
1, 2, 3, 4, 6, 7, 8, 9, 10 [polynomials].
basis functions
SINDy (esc)
[polynomials, sin, cos, exp]
optimizer threshold
0.05, 0. 0.15 0.025, 0.05, 0.075 20, 100
optimizer alpha optimizer max iterations
polynomial degree optimizer threshold
2, 3, 4, 5, 6, 7, 8, 9, 10
0.05, 0.1, 0.15
SINDy (poly)
optimizer alpha optimizer max iterations
0.025, 0.05, 0.075
20, 100
23

---

Published as a conference paper at ICLR 2024
1.0
1.0
1.0
Accuracy
Accuracy
Accuracy
0.5
0.5
0.5
odeformer (opt)
odeformer
0.0
10
20
30
10
20
30
10
20
30
T1 To
T1 - To
T1 F1-To
(a) p 0,
(b) 0, 0.01
(c) p 0.1,
0.01
1.0
1.0
1.0
median R
median
median
0.5
0.5
0.5
0.0
0.0
0.0
20
30
10
20
30
10
20
30
To
Ti To
T1-To
(d) p = 0
(e) p = 0, 0.01
(0 0.1, 0.01
Figure 8: Evaluation on Strogatz with missing chunks of data. We compare ODEFormer with and without constant optimization (opt). Originally, trajectories in the Strogatz dataset comprise 100 time points. For this evaluation drop all data in the interval [Tu Ti. and additionally drop fraction of p of the remaining time points uniformly at random. For σ > 0 the trajectories are additionally contaminated with noise. Plots in the upper row show results in terms of accuracy (fraction of samples with n2 > 0.9) whereas the lower row shows the median R22 score across the dataset.
F INTERPOLATION EXPERIMENTS FOR MISSING DATA CHUNKS
Whereas in the main text considered irregular grids and corruptions by Gaussian noise, here additionally evaluate ODEFormer in case of missing chunks of data. This scenario may for instance correspond faulty sensor which stops recording data for duration of time, e.g. until gets fixed or replaced.
To simulate this scenario use the Strogatz dataset. Each trajectory in this dataset covers the time interval [0, 10] and is originally comprised of 100 time points. Out of these 100 time points, we drop all time points within the interval [To, Ti] for [To, Ti] {[45, 55), [40, 60). [35, 65]). i.e., we assess three progressively more challenging missing data scenarios. Additionally, we corrupt the remaining data by randomly dropping observation and noise contamination as described in Section 5. We evaluate the same (trained) model as in Section 5. We remark that this model has not been optimized for this particular scenario as during training time points were dropped uniformly at random across the entire trajectory. Thus, intervals of entirely missing data may have been encountered in during training in exceptionally rare cases only.
model performance
Results under different signal corruptions are presented in Figure 8 while exemplary trajectories and corresponding predictions are displayed in Figure 9. Even under this challenging scenario, ODEFormer still performs reasonably well for the majority of samples as exemplified by both accuracy and median R² scores.
24

---

Published as a conference paper at ICLR 2024
R 0.986
R = 0.958
55
55
state
state
5
2
8
10
time
time
0.977
0.990
state
state
4
10
time
R4 0.850
R 0.959
18
state
state
0
2
8
10
0
2
4
6
10
time
time
R = 0,490
R = 0.987
9
state
state
-1
-1
0
4
8
10
2
6
10
time
time
R 0.467
R 0. ..594
12
12
state
state
-1
1
10
time
time
R 0.843
0.979
state
state
-3
-3
0
8
10
2
6
10
time
time
R< 0.809
R
0.843
2
state
state
1 0
2
4
6
8
10
0
2
6
8
10
time
time
Figure 9: Predictions of ODEFormer (left column) and ODEFormer (opt) (right column) noisy observations (black markers) when data is missing in a consecutive interval (red area). The thin gray line corresponds to the true (but unobserved) solution trajectory, the red line corresponds the trajectory given by integrating the predicted ODE, Each row corresponds to one of the seven ODEs available in dataset Strogatz.
25

---

Published as a conference paper at ICLR 2024
G EFFECT OF BEAM SIZE
Figure 10, we study the impact of the beam size on reconstruction and generalization performance. While reconstruction improves with the beam size, generalization hardly changes. This highlights the importance of using both metrics: the two are not necessarily correlated, and the latter is a much better proxy of symbolic recovery than the former.
1.00
Figure 10: Increasing the beam size im- proves reconstruction, but not generaliza- tion. We plot the average reconstruction and generalization RR-score 10,000 noise-free free, densely samples synthetic examples for vari- ous beam sizes and temperature of 0.1.
fitting
R2 score
0.75
generalization beam size 1 beam size
0.50
0.00
4 dimensions
H DETAILED COMPARISON WITH CLOSELY RELATED MODELS
H.1 SIMILARITIES AND DIFFERENCES
As mentioned Section 2, ODEFormer is conceptually similar to the models proposed by A'ooliii et al. (2022) and Becker et al. (2023). Like ODEFormer, these models leverage large-scale training on synthetically generated samples to learn a transformer model that maps numerical input sequences to symbolic output sequences. While the output for ODEFormer and NSODE (the model presented by Becker et al. (2023)) correspond to (the right-hand side of) symbolic ordinary differential equations (ODEs), the model by d' Ascoli et al. (2022) predicts the discrete recurrence relation of the input sequence. A major difference between ODEFormer and these two related models is that ODEFormer works with multivariate inputs, i.e., it is able handle systems of differential equations, whereas both other models only support univariate input sequences. For this reason neither NSODE nor the model by d' Ascoli et al. (2022) include an embedding module that needs combine multivariate inputs as described in Section 4.
While the recurrent sequence prediction approach outlined by d' A'scolii et al. (2022) already delineates a distinct use-case, the divergence between ODEFormer and NSODE is especially noteworthy from a practical perspective. While multivariate systems of ODEs can capture rich dynamic behaviors (including e.g. limit cycles), univariate ODEs are somewhat limited in that they can only capture monotonically converging or monotonically diverging behavior. Moreover, the majority of intricate modeling scenarios typically involves multiple dynamically interacting variables, rendering univari- ate model simply inadequate. Consequently, ODEFormer significantly broadens the practical toolset for ODE modelers in practice. We list the technical differences between all three models in Table 6.
Table 6: Model architecture of ODEFormer in comparison to closely related models.
ODEFormer (ours)
Becker et al. (2023)
d'Ascoli et al. (2022)
multivariate ODE
discrete univariate
Developed for
univariate ODEs
systems
recurrence sequences
layers
Encoder
attn. heads
16
16
8
embedding dim. forward dim.
512
512
512
512
2048
512
positional enc. layers
learned
16
6
Decoder
# attn. heads
16
16
8
embedding dim.
512
512
512
forward dim. beam size constants enc.
512
2048
512
50
1536
sign, mantissa, exp.
two-hot
sign, mantissa, exp.
26

---

Published as a conference paper at ICLR 2024
H.2 PERFORMANCE COMPARISON
Experimental setup. As NSODE can not be applied to multivariate systems, we can not evaluate its performance on Strogatz and ODEBench. Instead, we evaluate ODEFormer the datasets used in Becker et al. (2023), following their exact experimental conditions, which we briefly restate here for self-containment:
The evaluation includes three datasets: Textbook contains 12 ODEs transcribed from physics textbooks, Classic contains 26 ODEs which correspond classic (algebraic) equations in the symbolic regression literature which are re- re-iteerreeted as ODEs and lastly Large comprises 163 ODEs which are synthetically generated following the algorithm of Lample & Charton (2019). Although this dataset is generated following similar procedure described in Section 3, we expect the distribution of the generated data be quite different from ODEFormer's training distribution essentially all syper-parameters (c.g. maximum. tree depth, operator probabilities, probability to sample variable vs. a constant) differ. For details on the construction of the Large dataset, section 3.1 in Becker et al. (2023). All equations are integrated scipy. integrate odeint (Virtanen et al., 2020) to interface the LSODA (Hind- marsh & Laboratory, 1982) solver package. Numerical solution are subsequently systematically randomly subsampled, resulting in an equidistant irregular sampling grid. We remark that trajectories with 256 time points contain more time points than ODEFormer has seen for any trajectory during training (max. 200, see Section 3). Furthermore the numerical solver differs from the one used for generating the training data according to Section 3. Trajectories are corrupted with multiplicative noise drawn from A(1,
the interval [0, 2] with 512 steps using
Evaluation always occurs on the original non-subsampled and noise-free free solution trajectory.
Results for ODEFormer are obtained with beam size of 50 and without additional constant opti- mization. Performance is evaluated in terms of median Rt (which is the main metric used in Becker et al. (2023)) as well as accuracy as defined in Section 5, that is, fraction of predictions for which R 0.9.
The results in Tables 7 to 9 show that across all three datasets and all experimental conditions, ODEFormer performs quite similar NSODE. In particular on Textbook and Large both models achieve strong results with negligible differences in median R²2 yet slightly superior performance of NSODE terms of accuracy. This trend is also observable on dataset Large where differences in median R are again small in most cases yet differences in accuracy are more pronounced favor NSODE. This result expected given that equations in Large are generated from the same distribution that NSODE is trained on while the training distribution for ODEFormer differs. Overall, we conclude that ODEFormer generalizes well to different test time distributions.
Results.
Table 7: Results on Textbook
model
#points
sampling
a (noise)
median R²
accuracy
ODEFormer
128
equidistant equidistant
0.99
NSODE
128
0
0.99
ODEFormer
128
irregular
0.01
0.99
0.91
NSODE
128
irregular
0.01
0.999
ODEFormer
128
inregular
0.02
0.994
0.818
NSODE
128
inregular
0.02
0.988
0.91
ODEFormer
256
equidistant
0
0.999
NSODE
256
equidistant
O
0.997
ODEFormer
256
inregular
0.01
0,994
NSODE
256
irregular
0.01
0.999
ODEFormer NSODE
256
irregular
0.02
0.993
256
irregular
0.02
0.987
27

---

Published as a conference paper at ICLR 2024
Table 8: Results on Classic
model
#points
sampling
(noise)
median R²
accuracy
ODEFormer
128
equidistant
0.999
0.92
NSODE
128
equidistant
0
0.999
ODEFormer
128
irregular
0.01
0.998
0.92
NSODE
128
irregular
0,01
0.999
ODEFormer
128
irregular
0.02
0.998
0.84
NSODE
128
inregular
0.02
0.997
0.88
ODEFormer
256
equidistant
0
0,999
0.92
NSODE
256
equidistant
O
0.999
ODEFormer
256
inregular
0,01
0,998
0,88
NSODE
256
irregular
0.01
0.999
ODFForme NSODE
256
irregular
0.02
0,998
0.84
256
irregular
0.02
0.993
0.92
Table Results on Large
model
#points
sampling a (noise)
median R²
accuracy
ODEFormer
128
equidistant equidistant
0.978
0.652
NSODE
128
0
0.998
0.789
ODEFormer
128
irregular
0.01
0.961
0.571
NSODE
128
irregular
0.01
0.997
0.776
ODEFormer
128
inregular
0.02
0,933
0.565
NSODE
128
irregular
0.02
0.982
0.671
ODEFormer
256
equidistant
0
0,996
0.64
NSODE
256
equidistant
0
0.998
0.801
ODEFormer
256
irregular
0,01
0.99
0.627
NSODE
256
irregular
0.01
0.999
0.82
ODEFormer
256
irregular
0.02
0.923
0.522
NSODE
256
irregular
0.02
0.991
0.745
28

---

Published as a conference paper at ICLR 2024
I ADDITIONAL RESULTS ON BENCHMARKS
In Figures 11 to 16, we plot histograms to better visualize how each model performs across different noise conditions, datasets, and evaluation tasks.
ODEFormer topt)
ODEFormer (opt)
ODEFormer (opt)
ODEFormer
20
ODEFormer
20 GDEFormer
0
PySR
PySR
PYSR
SINDY (pdly)
20
SINDY (poly)
10
SINDY (paly)
20
10
ENC
ENC
EHC
SINDY (eac)
20
SINDY (esc)
20
SINDy (eec)
0-
5-
(pol))
ProGED (pol))
10
ipoly)
0
PreGED
ProGED
ProGED
20
SINDY
28 SINDY
28 SINDY
0
EPLEX
20
EPPLX
28
EPLEX
20
28
0
Inalld
<0
0.0
0.2
0.4
0.6
0.8
1.0
Inwild
<0
no
0.2
044
0.6
0.8
1.0
Inwali
0.0
0.2
0.4
0.6
0.8
T.
8²
R²
R?
(a) o=0
(b) 0=0.01
a=0.05
Figure 11: Histogram of per equation R² scores for the reconstruction task on Strogatz. Subfigures correspond to different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions as well as the number predictions that yielded an It2 score below 0. The red dashed line corresponds to the mean R² score across equations, the black dashed line corresponds to the median R² score.
29

---

Published as a conference paper at ICLR 2024
20
ODEFormer (opt)
20
(opt)
20
ODEFormer (opt)
0
ODEFormer
ODEFormer
ODEFormer
PYSR
20 PySR
10
PySR
20
SINDY (poly)
20 SINDY (poly)
20
SINDy
0
FFX
FFX
20
EHC
20
EHC
10
EHC
0
SINDY
SINDY
20
SINDy
PreGeo (poly)
ProteDD (pdly)
ProGED (paly)
ProGen
20
ProGED
20
ProGED
0
SINDY
28
AFR
20
FE AFP
20
FE ARP
20
28
0
AFP
Inalid
0.0
0.2
0.5
0.9
1.0
Involld
0.0
0.2
0.6
0.8
1.0
Iwalld
0.0
0.2
0.4
0.8
1.0
(a) o=0
(b) 0=0.01
(c) a=0.05
Figure 12: Histogram of per equation R² scores for the reconstruction task on Strogatz where 50% % of the trajectory are dropped uniformly at random (p 0.5). Subfigures correspond to different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions as well as the number of predictions that yielded an It2 score below 0. The red dashed line corresponds the mean R score across equations, the black dashed line corresponds to the median R² score.
30

---

Published as a conference paper at ICLR 2024
ODEFormer
ODEFormer
ODEFormer
ODEFormer
PySR
63 PySR
40
PySR
0
FFX
SINDY (poly)
40
SINDY (poly)
20
sinty
ProGED (poly)
40
(pdly)
40
(paly)
0
EMC
ProGED
40
ProdED D
40
ProGEe
0
SINDY (esc)
SINDY (esc)
40
SIND
0
AFP
AFP
AFP
SIND
50 SINDY
63 SINDY
0
FE-ARR
EPLEX
EPLEX
EPLEX
Inwalid
Iwalld
(a) o=0
(b) 0=0.01
(c) 0=0.05
Figure 13: Histogram of per equation Ri2 scores for the reconstruction task on ODEBench. Subfig- ures correspond to different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions score below 0. The red dashed line corresponds the mean R2 score across equations, the black dashed line corresponds to the median R² score.
well as the number of predictions that yielded an It
31

---

Published as a conference paper at ICLR 2024
50
ODEFormer (opt)
50
(opt)
50
ODEFormer (opt)
0
ODEFormer
ODEFormer
ODEFormer
PYSR
40
PySR
40
PySR
20
40
0
SINDY
SINDY (poly)
SINDy
ProGED (poly)
40
ProGED (pol))
40
ProGED (poly)
0
ProGeD
ProGED D
ProGEe
50
SINDY (esc)
40
SINDY (esc)
63
SiND
0
SINDY
50
SINDY
63
SINDY
40
40
0
EPLEX
EPLEX
Inalld
0.0
0.2
0.5
0.9
1.0
Involld
0.0
0.2
0.6
0.8
1.0
Iwaild
0.0
0.2
0.4
0.8
1.0
(a) a=0
(b) 0=0.01
(c) a=0.05
Figure 14: Histogram of per equation R²2 scores for the reconstruction task on ODEBench where 50% % of the trajectory are dropped uniformly at random (p 0.5). Subfigures correspond to different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions as well as the number of predictions that yielded an R² score below O. The red dashed line corresponds to the mean R² score across equations, the black dashed line corresponds the median R² score.
32

---

Published as a conference paper at ICLR 2024
ODEFormer (opt)
ODEFormer (opt)
ODEFormer (opt)
FFX
40
FFX
40
FFX
0
SINDY
SINDy
SINDy
PySR
40
PySR
40
PySR
(esc)
40
SINDgress)
63
SINDy
0
ProGep
ProGee
ProGED
EHC
50
EHC
40
EHC
0
50
0
FE-ARM
FE-AED
FE-AR
AFP
50
AFP
03: AFP
0-
SINDY
EPLEX
EPLEX
EPLEX
Inwilid
Iwalld
(a) o=0
(b) 0=0.01
(c) 0=0.05
Figure 15: Histogram of per equation R²2 scores for the generalization task on ODEBench. Subfig- ures correspond to different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions score below 0. The red dashed line corresponds the mean R2 score across equations, the black dashed line corresponds to the median R² score.
well as the number of predictions that yielded an It
33

---

Published as a conference paper at ICLR 2024
ODEFormer
40
40 ODEFormer
0
ODEFormer
ODEFormer
ODEFormer
FFX
40
FFX
50
FFX
SINDy
40
SINDy
50
SIND
0
PYSR
PySR
PySR
40
(esD
50
03
SINDw
0
(poly)
(pdl))
50
Ipoly)
EHC
EHC
EHC
50
ProGED
50
ProGe
50 ProGen
0
FE-AFD
AFP
50
AFP
63
AFP
50
63
SINDY
63
SINDY
0
EPLEN
EPLEX
EPLEXI
Inalld
0.0
0.2
0.5
0.9
1.0
Involld
0.0
0.2
0.6
0.8
1.0
Iwalld
0.0
0.2
0.4
0.8
1.0
(a) a=0
(b) a=0.01
(c) a=0.05
Figure 16: Histogram of per equation 2²2 scores for the generalization task on ODEBench where 50% of the trajectory different noise levels. The y-axis represents counts and is scaled per model for better visibility of the distribution of scores. The x-axis annotations "invalid" and "<0" respectively denote the number of invalid predictions as well as the number of predictions that yielded an R² score below O. The red dashed line corresponds to the mean R² score across equations, the black dashed line corresponds the median R² score.
dropped uniformly at random (p 0.5). Subfigures correspond to
34


Reviews:
[{"review_id": "HSWv9egCvL", "reviewer": "Reviewer_hbTD", "summary": "The paper studies the problem of symbolic regression for dynamical systems, specifically ODE, with the use of transformers. Apart from the adjustments need to use a transformer based model for this task, authors also propose a new benchmark. This is claimed to be more diverse and larger than existing ones. Empirically, the proposed method is shown to outperform existing baseline methods in terms of reconstruction as well as generalization.", "questions": "1) Please compare and contrast the proposed method with the work \"Becker, Sören, et al. \"Predicting Ordinary Differential Equations with Transformers.\" (2023).\". As acknowledged that this is closely related, it should be clearly articulated what are the differences. The mentioning of difference between univariate and multivariate is okay, but this needs more explanation. I strongly believe this should be a baseline for comparison in some setting.\n\n2) How should one think about inferring PDE with a variant of this framework, can we have a discussion around this as part of limitation if that is the case? \n\n3) The learning from multiple trajectories, doesn't seem to give promising performance, do authors have comments on explaining this behavior?\n\n4) For the noise corruption of the data, the authors have tried out adding Gaussian noise and dropping samples uniformly at random. Can authors comment on missing chunks of data instead. This might be of practical implication in certain cases, where the data collection mechanism (sensor) went faulty for some time before resuming activity.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is extremely well-written (including notations, clearly stating contributions, etc.), very well motivated and the proposed method is shown to achieve state of the art performance.\n\nThe placement within existing literature is very well articulated.\n\nSince, authors propose a benchmark, it is appreciated that the data generation procedure is outlined precisely.\n\nThe section on filtering of the data to avoid rapidly converging and divergent systems is worth mentioning, details like these make a benchmark standout.\n\nTokenization, embedding process and encoding of the symbolic functions is very well justified.\n\nBaseline methods have been chosen appropriately and the empirical evidence is very convincing.", "weaknesses": "The authors have mentioned a few of their limitations, which is great. \n\nWhile authors mention the presence of a very related work \"Becker, Sören, et al. \"Predicting Ordinary Differential Equations with Transformers.\" (2023).\", I am not sure why they do not elucidate the difference between the paper mentioned and their proposed method, why is this method not used as a baseline? There is no attempt to illustrate the difference at all.\n\nI have another point which I would like to bring up with the authors regarding the generation of the data, specifically the way it is being integrated. As mentioned (and I think this is reasonable as a first step) the authors using fixed homogeneous grid to integrate, and it is claimed that number of points don't matter during inference (Figure 3), is this indeed an artifact of the way the data is generated. How do the results change when the integration procedure is altered, it would be great to have the authors comment on this issue.\n\nPlease see questions section for further.", "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": "", "extra_scores": {}}, {"review_id": "MQJJiHbUKc", "reviewer": "Reviewer_znUh", "summary": "ODEFormer, a transformer model designed for dynamical symbolic regression. This model is capable of inferring multidimensional ordinary differential equation (ODE) systems from noisy and irregularly sampled data. It utilizes a pre-trained sequence-to-sequence transformer on synthetic data to generate symbolic expressions directly from observations. The paper also introduces ODEBench, a benchmark dataset for dynamical symbolic regression, comprising 63 ODEs sourced from the literature, modeling real-world phenomena across dimensions one to four, including chaotic systems.\n\nThe reviewer further highlights the evaluation and comparison of ODEFormer with existing methods. It assesses ODEFormer's performance on both the established Strogatz dataset and the new ODEBench dataset. The comparison encompasses various techniques based on genetic programming, regression, and Monte Carlo methods.", "questions": "One key aspect that warrants further exploration is the practical application of the technique on real-world datasets. While the framework shows promise in artificially constructed datasets, its effectiveness in solving real-world problems, where data can be noisy, irregularly sampled, and complex, remains to be demonstrated. Evaluating its performance on non-artificial, real-world datasets across various domains, such as finance, healthcare, or environmental monitoring, would provide valuable insights into its applicability and limitations in practical scenarios.\n\nTime series data indeed represents a compelling use case for the framework. Time series forecasting is a critical application in various fields, including finance, energy, and climate modeling. Therefore, it is crucial to investigate the applicability of this method in a time series forecasting modeling scenario. Demonstrating its effectiveness in accurately modeling and predicting time-dependent data can significantly enhance its practical utility and establish its relevance in solving real-world, dynamic data challenges.", "limitations": "", "rating": 3, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 1, "strengths": "The introduction of ODEFormer represents a novel framework that serves the purpose of generating ordinary differential equations (ODEs) specifically designed for testing dynamical systems. This innovative approach allows researchers and practitioners to create ODE models that can accurately capture and represent the dynamics of real-world systems, providing a valuable tool for testing and understanding the behavior of complex dynamical systems.\n\nODEFormer introduces a pioneering use case for transformers in the realm of ODEs. Traditionally, transformers are employed in natural language processing and sequential data tasks. However, in this context, they are harnessed to directly infer ODE systems from noisy and irregularly sampled data. This expansion of transformer applications into the domain of ODEs signifies a breakthrough, offering a versatile and data-driven approach for modeling and analyzing complex dynamic systems, and opening up new possibilities for the fusion of machine learning techniques with physics-based modeling.", "weaknesses": "One notable limitation in the presented work is the absence of comparisons with benchmarks employed in previous studies, such as the widely recognized benchmark datasets used in Neural ODE (Chen et al). This absence makes it challenging to gauge how the proposed ODEFormer framework performs in comparison to existing approaches on well-established and widely accepted testing scenarios.\n\nThe demonstrated applicability of ODEFormer on toy datasets represents another potential limitation. Toy datasets are typically simplistic and may not fully capture the complexity and variability encountered in real-world applications. Therefore, the extent to which ODEFormer can effectively handle and model more complex, real-world data remains an open question and warrants further investigation and validation on diverse and challenging datasets.", "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": "", "extra_scores": {}}, {"review_id": "ulp7r9iBT3", "reviewer": "Reviewer_hEs1", "summary": "The authors proposed a model for discovering the governing law of dynamical systems out of observed trajectory. The model treats trajectories as sequences of tokenized numbers and utilizes a transformer network to learn the symbolic form of governing laws from them. The authors conduct the experiment on 1D to 4D systems, and the proposed method outperforms the baseline methods.", "questions": "(1)\tCould the authors give a detailed description and specific form of the loss function used during training? Is it symbolic or numerical? If the loss function is symbolic, does it measure the similarity between the symbolic tree structure or simply between the sequences? And as the ground-truth equations can be organized in different orders under both forms, which one of them is chosen in practice, and what are the reasons for refusing other possible sequences/trees?\n\n(2)\tAs the authors claimed in section 2, most existing approaches of symbolic regression require a separate optimization for each new observed system. However, it is unclear in the paper whether or not the proposed model needs separate optimization on unseen systems as previous models do. Could you add some analysis on this problem and provide reasons and results to support your arguments?\n\n(3)\tThe authors conduct subsampling in the experiment part and compare model performances under different subsampling rates. Could you explain further what condition in practice the subsampling process corresponds to? And how does the proposed model deal with the resulting irregular time-interval as most sequence transformers are built to encode trajectories with regular time intervals?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "(1) The paper is well motivated, as discovering the symbolic form of governing laws from observed data has always been the focus of scientific research.\n\n(2) The paper contributes to neural network models for symbolic regression of ordinary differential equations.", "weaknesses": "(1)\tThe idea in this paper is not novel enough. The idea of transformer-based symbolic regression on tokenized sequences was previously reported[1]. The proposed model shares the same network structure but is applied to ordinary differential equation data. However, no analysis could be found in the paper on how the network is adapted to this new type of data.\n\n(2)\tFigures and explanations of the proposed model are way too rough and lack the necessary details.\n\n(3)\tThis paper lacks crucial detail on dataset separation and the definition of tasks in the experiment section.\n\nReference:\n\n[1] d'Ascoli S, Kamienny P A, Lample G, et al. Deep symbolic regression for recurrent sequences[J]. arXiv preprint arXiv:2201.04600, 2022.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)