End-to-end Symbolic Regression with Transformers
Authors: Pierre-Alexandre Kamienny
Stéphane d'Ascoli
Guillaume Lample
Francois Charton
Venue: NeurIPS 2022 Accept (2022)
Keywords: symbolic regression
transformers
supervised learning
deep learning

Abstract: Symbolic regression, the task of predicting the mathematical expression of a function from the observation of its values, is a difficult task which usually involves a two-step procedure: predicting the "skeleton" of the expression up to the choice of numerical constants, then fitting the constants by optimizing a non-convex loss function. The dominant approach is genetic programming, which evolves candidates by iterating this subroutine a large number of times. Neural networks have recently been tasked to predict the correct skeleton in a single try, but remain much less powerful.

In this paper, we challenge this two-step procedure, and task a Transformer to directly predict the full mathematical expression, constants included. One can subsequently refine the predicted constants by feeding them to the non-convex optimizer as an informed initialization. We present ablations to show that this end-to-end approach yields better results, sometimes even without the refinement step. We evaluate our model on problems from the SRBench benchmark and show that our model approaches the performance of state-of-the-art genetic programming with several orders of magnitude faster inference. 

End-to-end Symbolic Regression with Transformers
Pierre-Alexandre Kamienny" 1.2, Stéphane d' A'Acoll' 1,3
Guillaume Lample' François Charton
Meta AI
FSSR MIIA, Sorbonne Université
3 Department of Physics, Ecole Normale Supérieure pakami enwy@meta.ccccooommm
Abstract
Symbolic regression, the task of predicting the mathematical expression of a function from the observation of its values, is a difficult task which usually involves two-step procedure: predicting the "skeleton" of the expression up the choice of numerical constants, then fitting the constants by optimizing non-convex loss function. The dominant approach genetic programming, which evolves candidates by iterating this subroutine large number of times. Neural networks have recently been tasked to predict the correct skeleton a single try, but remain much less powerful.
In this paper, challenge this two-step procedure, and task a Transformer directly predict the full mathematical expression, constants included. One can subsequently refine the predicted constants by feeding them the non-convex optimizer as an informed initialization. We present ablations to show that this end-to-end approach yields better results, sometimes even without the refinement step. We evaluate our model problems from the SRBench benchmark and show that our model approaches the performance of sate-oofthh-an genetic programming with several orders of magnitude faster inference.
Introduction
Inferring mathematical laws from experimental data is central problem in natural science; having observed a variable y at n points [a,Jve,, implies finding a function such that y/s Rs f(x;,) for all i N,. Two types of approaches exist to solve this problem. In parametric statistics (PS), the function ff is defined by a small number of parameters that can directly be estimated from the data. On the other hand, machine learning (ML) techniques such decision trees and neural networks select f from large families of non-linear functions by minimizing a loss over the data. The latter relax the assumptions about the underlying law, but their solutions more difficult interpret, and tend to overfit small experimental data sets, yielding poor extrapolation performance.
Symbolic regression (SR) stands middle ground between PS and ML approaches: f is selected from a large family of functions, but required to be defined by an interpretable analytical expression. It has already proved extremely useful in variety of tasks such inferring physical laws 2
SR is usually performed in two steps. First, predicting a "skeleton". parametric function using pre-defined list of operators typically. the basic operations (+, x, +) and functions (sqrt, exp, sin, etc.). It determines the general shape of the law up a choice of constants, e.g. f(z cos(ar + b). Then, the constants in the skeleton (a, b) are estimated using optimization techniques, typically the Broyden Broyden-Fletcherr Goldfarb- Shanno algorithm (BFGS).
equal contribution
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Inference time (seconds)
Sca GP)
10
BRLE
10
05
AFF (37M4T))
10a
102
Ours
Ours (no ref)
Ours (skel)
10a
10
(skel)
Ours
Linn
27X
Ours tno reD
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Black-box median R
Feynman mean accuracy >R²
DL based SR
GPbaaed SR
ML methods
Figure Our model outperforms previous DL--based methods and offers least an order of magnitude inference speedup compared to SOTA GP-based methods. Pareto plot comparing the average test performance and inference time our models with baselines provided by the SRbench benchmark [7]. both Feynman SR problems [ and black-box regression problems. We use colors to distinguish three families of models: deep-learning based SR, genetic programming- -aaedd SR and classic machine learning methods (which do not provide symbolic solutions). A similar Pareto plot against formula complexity is provided in Fig.
The leading algorithms for SR rely genetic programming (GP). At each generation, a population of candidates is predicted, and the fittest ones are selected based on the data, and mutated to build the next generation. The algorithm iterates this procedure until satisfactory accuracy level is achieved.
While GP algorithms achieve good prediction accuracy, they are notably slow (see the Pareto plot of Fig. D. Indeed, the manually predefined function space to search is generally vast, and each generation involves costly call to the BFGS routine. Also, GP does not leverage past experience: every new problem is learned from scratch. Inference time, i.e. time required to output satisfactory expression, for most GP algorithms is both long and unbounded (the longer, the better the results), therefore this property has been neglected by the SR community, however fast inference can be useful to improve search-basee based SR algorithms, e.g. by reducing search space or providing initial guesses, as well as tackle practical applications with time constraints, e.g. control reinforcement learning 3
Supervised training neural networks built for language modelling on large datasets of synthetic examples has recently been proposed for SR 5 (). These references follow the two-step procedure (predicting the skeleton then fitting the constants) inherited from GP. Once the model is trained, at inference, the skeleton is predicted via a simple forward pass, and a single call to BFGS is needed, thus resulting in significant speed-up compared to GP. However, these methods not accurate as state-of of-the-art GP, and have so far been limited to dideensiiiinal functions (D < 3). We argue that two reasons underlie their shortcomings.
First, skeleton prediction an ill-posed problem that does not provide sufficient supervision: different instances of the same skeleton can have very different shapes, and instances of very different skeletons can be very close. Second, the loss function minimized by BFGS can be highly non-nonconvex: even when the skeleton is perfectly predicted, the correct constants are not guaranteed to be found. For these reasons, we believe, and will show, that doing away with skeleton estimation as a intermediary step can greatly facilitate the task of SR for language models.
Contributions In this paper, we train Transformers over synthetic datasets to perform end-to-end (E2E) symbolic regression: solutions are predicted directly, without resorting skeletons. To this effect, we leverage a hybrid symbolie-numeric vocabulary, that uses both symbolic tokens for the operators and variables and numeric tokens for the constants. One can then perform a refinement of the predicted constants by feeding them as informed guess to BFGS, mitigating non- linear optimization issues. Finally, we introduce generation and inference techniques that allow models to scale to larger problems: up to 10 input features against 3 in concurrent works.

---

Evaluated over the SRBench benchmark [7]. our model significantly narrows the accuracy gap with state-of- -hetttr-tt GP techniques, while providing several orders of magnitude of inference time speedup (see Fig. D.. We also demonstrate strong robustness to noise and extrapolation capabilities.
Related work
SR challenging task that traces back from few decades ago, with a large number of open-source and commercial softwares, and has already been used to accelerate scientific I0 Most popular frameworks for symbolic regression use GP [ 12 [ 14
discoveries g
15
19 (see [7] for a recent review), but SR has also seen growing interest from the Deep Learning (DL. community, motivated by the fact that neural networks are good at identifying
16
17
18
qualitative patterns.
Neural networks have been combined with GP algorithms, e.g. to simplify the original dataset M. or to propose a good starting distribution over mathematical expressions[20) [21 22] propose modifications feed-forward networks to include interpretable components, i.c. replacing usual activation functions by operators such as cos, sin, however these are hard to optimize and prone to numerical issues.
Language models, and especially Transformers 23 have been trained over synthetic datasets to solve
various mathematical problems: integration [24] dynamical systems (25), linear algebra [26). formal
logic [27] and theorem proving 28 A few papers apply these techniques to symbolic regression: the aforementioned references 6 train Transformers predict function skeletons, while (2) infers
one-dimensional recurrence relations in sequences of numbers. (30] trains fully-connected networks to predict simple formulas from tabular data.
The recently introduced SRBench 7 provides a benchmark for rigorous evaluation of SR methods, in addition 14 SR methods and 7 ML baselines which will compare to in this work.
1 Data generation
Our approach consists in training language models on vast synthetic datasets. Each training example is a pair: a of N points (z, y) E R R as the input, and function f such that y = f(z) the targer Examples are generated by first sampling a random function /. then set of N input values (aN in D and computing y4 (1:)).
1.1 Generating functions
To sample functions f. we follow the seminal approach of Lample and Charton [24, and generate random trees with mathematical operators as internal nodes and variables or constants as leaves. The procedure is detailed below (see Table in the Appendix for the values of parameters):
1. Sample the desired input dimension D of the function f from (111, )
2. Sample the number of binary operators b from U&D - D + ) then sample b operators
from 24 (++,,,+
3. Build a binary tree with those nodes, using the sampling procedure of 24
4. For each leaf in the tree, sample one of the variables Fd, d E ND.
5. Sample the number of unary operators from (0, Umax then sample u operators from the list O2, in Table 3 and insert them random positions in the tree.
6. For each variable Ta and unary operator u, apply a random affine transformation, i.e. replace
Za by ar, + b, and u by au + with (a, b) sampled from Daft-
Note that since we require independent control on the number of unary operators (which is indepen- dent of D) and binary operators (which depends on D), we cannot directly sample a unary-binary tree as in (24]. Note also that the first D variables are sampled in ascending order to obtain the desired input dimension, which means functions with missing variables such as I + are never encountered; this is not an issue as our model can always set the prefactor of Tz to zero. As discussed quantitatively
We only consider functions from Ie into E; the general case f RR independent subproblems. "Note that although the division operation is technically binary operator, it appears much less frequently than additions and multiplications in typical expressions [3I]. hence we replace it by the unary operator
1 can be handled as P
inv: + 1/z.

---

Target
Inpit
ON. 2i11)
Cromentropy
Oulput
Embodder
Transormer
Scaled
Scaled
Training
Inference
Figure 2: Sketch of our model. During training, the inputs are all whitened. At inference, whiten them as a pre-processing step; the predicted function must then be unscaled to account for the whitening.
in App. C the number of possible skeletons as well as the random sampling of numerical constants guarantees that our model almost never sees the same function twice, and cannot simply perform memorization. See App. for examples of the skeleton of generated expressions.
1.2 Generating inputs
For each function f R. - R, we sample N E UM10D, N input values I, E R from the distribution Dr described below, and compute the corresponding output values y/ f(z)). If any I; is outside the domain definition of f or any y/ is larger 10a00, the process is aborted, and we function. Note that rejecting and resampling out-of-domain values of Ti the obvious and cheaper alternative, would provide the model with additional information about f. by allowing it to learn its domain of definition.
start again by generating
To maximize the diversity input distributions seen during training, we sample our inputs from a mixture of distributions (uniform or gaussian), centered around & random centroids fl see App. for some illustrations
Input samples are generated as follows:
1. Sample a number of clusters 1 ~ ((1,kmaa and k weights w; ~ U(0 1), which are then normalized so that Σ, w; =
2. For each cluster
Nx. sample a centroid ~ (0, a vector of variances o, ~ 24(0, 1) and distribution shape (gaussian uniform) D; E (N,U).
3. For each cluster E Nr. sample [o,NN] input points from Di(44, o;) then apply a random rotation sampled from the Haar distribution.
4. Finally, concatenate all the points obtained and whiten them by substracting the mean and dividing by the standard deviation along each dimension.
1.3 Tokenization
we represent numbers in base 10 floating-point notation, round them to four significant digits, and encode them as sequences of tokens: their sign, mantissa (between 0 and 9999), and exponent (from E-100 to E100).
Following 26]
To represent mathematical functions sequences, we enumerate the trees in prefix order, i.e. direct Polish notation, as in 24 operators and variables and integers are represented as single autonomous tokens, and constants encoded explained above.
For example, the expression fz) cos(2.4242.) encoded as [cos, muul, 2424,,-33xx0 Note that the vocabulary of the decoder contains mix of symbolic tokens (operators and variables) and numeric tokens, whereas that of the encoder contains only numeric token
'For + 0o, such a mixture could in principe approximate any input distribution. 'The embeddings of numeric tokens are shared between the encoder and decoder.

---

Head
Head 2
Head 3
Head 4
Head
Head 6
Head 7
Head 8
Layor
w/
w/
Layer
Figure 3: Attention heads reveal intricate mathematical analysis. We considered the expression f(z) = sin(r)/ /z, with N = 100 imput points sampled between -20 and 20 (red dots; the y-axis is arbitrary). We plotted the attention maps of a few heads of the encoder, which are x N matrices where the element (i represents the attention between point and point j. Notice that heads 2, and of the second layer analyze the periodicity of the function Fourier-like manner.
Methods
Below describe approach for end-to-end symbolic regression; please refer to Fig. f an illustration.
2.1 Model
Embedder Our model is provided N input points (x. E DR++1, each of which is represented as 3(D + 1) tokens of dimension demb- As and become large, this results in long input sequences (e.g. 6600 tokens for D 10 and N 200), which challenge the quadratic complexity of Transformers. To mitigate this, we introduce an embedder to map each input point to a single embedding.
The embedder pads the empty input dimensions Dmax then feeds the 3(Dmax + dedemb L-)ddmm-dimenstee dimensional vector into 2-layer fully-connected feedforward network (FFN) with ReLU activations, which projects down dimension deme 6 The resulting N embeddings of dimension demb are then fed the Transformer.
Transformer We sequence sequence Transformer architecture 23 with 16 attention heads and an embedding dimension of 512, containing a total of 86M parameters. Like [26]. we observe that the best architecture for this problem is asymmetric. with a deeper decoder: we use layers in the encoder and 16 in the decoder. A notable property of this task is the permutation invariance of the N input points. To account for this invariance, we remove positional embeddings from the encoder.
As shown in Fig. B and detailed in App. D the encoder captures the most distinctive features of the functions considered, such critical points and periodicity, and blends mix of short-ranged heads focusing on local details with long-ranged heads which capture the global shape of the function.
Training We optimize cross-entropy loss with the Adam optimizer, warming up the learning rate from 10-7 7 to 2.10 over the first 10,000 steps, then decaying it as the inverse square root of the number of steps, following (23) We hold out validation set of 10 examples from the same generator, and train our models until the accuracy on the validation set saturates (around 50 epochs of 3M examples).
Input sequence lengths vary significantly with the number points N: to avoid wasteful padding, we batch together examples of similar lengths, ensuring that full batch contains minimum of 10,000 tokens. 32 GPU with 32GB memory each, one epoch processed in about half an hour.
2.2 Inference tricks
In this section, we describe three tricks to improve the performance of our model at inference.
'We explored various architectures for the embedder, but did not obtain any improvement; this does not appear to be critical part of the model.
5

---

Table 1: The importance of an end-to-end model with refinement.
Model
Function f(z y)
Target
sin 10x)) exp(0.1y)
Skeleton + BFGS
(1. + 0.19)
E2E no BFGS
sin(9.92) exp(0.ly)
E2E + BFGS random init E2E BFGS model init
(0.0099)5 (0.20))7y)
sin(10x) cxp((.1y)
The skeleton approach recovers an incorrect skeleton. The E2E approach predicts the right skeleton. Refinement worsens original prediction when randomly initialized, and yields the correct result when initialized with predicted constants.
Refinement Previous language models for SR, such as [6], follow a skeleton approach: they first predict equation skeletons, then fit the constants with non-linear optimisation solver such as BFGS. In this paper, we follow end-to-end (E2E) approach: predicting simultaneously the function and the values of the constants. However, we improve our results by adding a refinement step: fine-tuning the constants posteriori with BFGS, initialized with our model predictions"
This results in large improvement over the skeleton approach, as we show by training Transformer to predict skeletons in the same experimental setting. The improvement comes from two reasons: prediction of the full formula provides better supervision, and helps the model predict the skeleton; second, the BFGS routine strongly benefits from the informed initial guess, which helps the model predict the constants. This is illustrated qualitatively in Table and quantitatively in Table g
first,
Scaling As described in Section 1.2. all input points presented to the model during training are whitened: their distribution is centered around the origin and has unit variance. To allow accurate prediction for input points with a different mean and variance, we introduce a scaling procedure during inference. Let f the function be inferred, x be the input points, and μ mean(x), a std(x). As illustrated in Fig. 2 we pre-process the input data by replacing æ by -_ The model then predicts f(i) = f(r + μ). and we can recover an approximation of f by unscaling the variables in f.
a
This gives our model the desirable property to be insensitive to the scale of the input points: DL-based approaches to SR are known to fail when the inputs are outside the range of values seen during training (29 26]. Note that here, the scale of the inputs translates to the scale of the constants in the function f: although these coefficients are sampled in Dan during training, coefficients outside Dair can be expressed by multiplication of constants in Dar-
Bagging and decoding Since our model was trained on N < 200 input points, it does not perform satisfactorily inference when presented with more than 200 input points. To take advantage of large datasets while accommodating memory constraints, we perform bagging: whenever N is larger than 200 inference, we randomly split the dataset into B bags of 200 input point:
For each bag, we apply a forward pass and generate function candidates via random sampling or beam search using the next token distribution. As shown in App. F [6],, the more commonly used beam search (34] strategy leads to much less good results than sampling due to the lack of diversity induced by constant prediction (typical beams will look like sin sin(1.1.) sin(0.9.),..... This provides us with a set of BC candidate solutions.
Inference time on one hand (which can be parallelized up to memory limits of the GPU), and the refinements of candidate functions on the other (which are CPU-based and could also be parallelized, although did not consider this option here).
Our model inference speed has two sources: the forward passes described above
"To avoid BFGS having to approximate gradients via finite differences, we provide the analytical expression the gradient using sympytorch (32) and functorch (3 Smarter splits, e.g.
could be envisioned, but not considered here.

---

Table 2: Our approach outperforms the skeleton approach.
Model
R²
Acco.1
Accp.01
AcCo.001
Skeleton + BFGS
0.43
0.40
0.27
0.17
E2E no BFGS
0.62
0.51
0.27
0.09
E2E + BFGS random init
0.44
0.44
0.30
0.19
E2E + BFGS model init
0.68
0.61
0.44
0.29
Metrics are computed over the 10, 000 examples of the evaluation set.
Since BC can become large, we rank candidate functions (according to their error on all input points). get rid of redundant skeleton functions and keep the best K candidates for the refinement step To speed the refinement, we use a subset of at most 1024 input points for the optimization. The parameters B, C' and K can be used as cursors in the speed-accuracy tradeoff:: in the experiments presented in Fig. we selected B 100, C = 10, K = 10.
3 Results
In this section, we present the results of our model. We begin by studying in-domain accuracy, then present results out-of-domain datasets.
3.1 In-domain performance
We report the in-domain performance of our models by evaluating them on fixed validation set of 100,000 examples, generated per Section [ Validation functions uniformly spread out over three difficulty factors: number of unary operators, binary operators, and input dimension. For each function, we evaluate the performance of the model when presented N [50, 100, 150, 200] input points (I, y), and prediction accuracy is evaluated on Mest 200 points sampled from fresh
instance of the multimodal distribution described in Section 1.2
We assess the performance of our model using two popular metries: R score 7 and accuracy to tolerance 6 29]
R==1 Σ Nea (yi )2)²
y1
Nua (y/ y)²
Acc,
max 15i5
5T)
Vi
where 1. is the indicator function.
R2t is classically used in statistics, but is unbounded, hence single bad prediction can cause the average Rt2 over set of examples to be extremely bad. To circumvent this, we set R22 = 0 upon pathological examples in (such examples occur in less that 1% of cases 10] The accuracy metric provides a better idea of the precision of the predicted expression as it depends on desired tolerance threshold. However, due to the presence of the max operator, it is sensitive outliers, and hence to the number of points considered at test time (more points entails a higher risk of outlier). To circumvent this, we discard the 5% worst predictions, following (6)
End-to-end outperforms skeleton In Table 2 we report the average in in-domain results of our models. Without refinement, our E2E model outperforms the skeleton model trained under the same protocol in terms of low precision prediction (Rt² and Acco.1 metrics), but small errors in prediction of the constants lead lower performance at high precision (Acca.001 metric). The refinement procedure alleviates this issue significantly, inducing three-fold increase in Acco.001 while also boosting other metrics. Initializing BFGS with the constants estimated in the E2E phase plays crucial role: with random initialization, the BFGS step actually degrades E2E performance. However, refinement with random initialization still achieves better results than the skeleton model: this suggests that the E2E model predicts skeletons better that the skeleton model.
are no guarantees that these would the after refinement, especially as optimization is particularly prone to spurious local optimas. Note that predicting the constant function y naturally yields H? score of O.
TThough these candidates are best functions without refinement,

---

1.0
1.0
A
B
C
0.8
0.8
Accuracy
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
0.0
2
10
Number of unary ops
Num of binary ops
Input dimension D
1.0
1.0
1.0
E
Interpolation
F
0.8
0.8
Extrapolation
0.8
Accuracy
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
0.0
50
100
150
200
10
10a
0.000
0.025 0.050 0.075 0.100
Number input pairs N
Test input variance
Noise variance
Skeleton
E2E ref
E2E random ref
E2E ref
Figure 4: Ablation over the function difficulty (top row) and input difficulty (bottom row). We plot the accuracy skeleton, E2E without refinement, E2E with refinement from random guess and E2E with refinement. A: number of unary operators. B: number of binary operators. C: input dimension. D: Low-resource performance, evaluated by varying the number of input points. E: Extrapolation performance, evaluated by varying the variance of the inputs. Robustness to noise, evaluated by varying the multiplicative noise added to the labels.
(Eq. D. see App. For for the R² score. We distinguish four models:
Ablation Fig. 4A,BCC presents an ablation over three indicators of formula difficulty (from left to right): number of unary operators. number of binary operators and input dimension. In all cases, increasing the factor of difficulty degrades performance, as one could expect. This may give the impression that our model does not scale well with the input dimension, but we show that our model scales in fact very well on out-of-domain datasets compared to concurrent methods (see Fig. 15 of the Appendix). We include a qualitative ablation on the improvement caused by the use of mixture of distributions in App.
Fig. 4D shows how performance depends the number of input points fed the model, N. In all cases, performance increases, but much more signicantly for the E2E models than for the skeleton model, demonstrating the importance of having a of data accurately predict the constants in the expression.
Extrapolation and robustness Fig. 4E, we examine the ability our models to interpo- late/cxtrapolate by varying the scale of the test points: instead of normalizing the test points to unit variance, we normalize them to a scale a. As expected, performance degrades as we increase a, however the extrapolation performance remains decent even very far away from the inputs (o 32).
Finally, in Fig. 4F, we examine the effect of corrupting the targets y with multiplicative noise of variance a: y + y(1 + E). ~ =). The results reveal something interesting: without refinement, the E2E model is not robust to noise, and actually performs worse than the skeleton model at high noise. This shows how sensitive the Transformer is to the inputs when predicting constants. Refinement improves robustness significantly, but the initialization of constants to estimated values has less impact, since the prediction of constants is corrupted by the noise.

---

3.2 Out-of-domain generalization
We evaluate our method on the recently released benchmark SRBench(7]. Its repository contains set of 252 regression datasets from the Penn Machine Learning Benchmark (PMLB) 335 in addition to 14 open-source and MI baselines. The datasets consist "ground-truth" problems where the true underlying function is known, as well as "black-box" problems which are more general regression datasets without an underlying ground truth. We filter out problems from SRBench to only keep regression problems with D <. 10 with continuous features; this results in 190 regression datasets, splitted into 57 black-box problems (combination of real-world and noisy, synthetic datasets), 119 SR datasets from the Feynman I and 14 SR datasets from the ODE-Strogatz [6] databases. Each dataset is split into 75% training data and 25% test data, on which performance is evaluated.
The overall performance of our models is illustrated in the Pareto plot of Fig. where we see that on both types of problems, our model achieves performance close to state-of-the-art GP models such as Operon with fraction of the inference time Impressively, our model outperforms all classic MI. methods (e.g. XGBoost and Random Forests) on real-world problems with lower inference time, and while outputting an interpretable formula.
We provide more detailed results on Feynman problems in Fig. 5 where we additionally plot the formula complexity, i.e. the number of nodes in the mathematical tree (see App. [ for similar results on black- box and Strogatz problems). Varying the noise applied the targets noise, we see that our model displays similar robustness to state-of-the-an GP models. We additionally include ablation on the use of scaling during inference in App. E
While the average accuracy or our model is only ranked fourth, outputs formulas with lower complexity than the top 2 models (Operon and SBP-GP), which is important criteria for SR problems: see App. TT for complexity-acouu aaccurrcy Pareto plots. To the best of our knowledge, our model is the first non-GP approach achieve such competitive results for SR.
Mean accuracy (R² 0.99)
Formula complexity
Inference time (seconds)
Operon
SBP-CP
GOMEM
Ours
EPLEX
MRCP
AFP FE
KG
Ours (skel)
AlFeynman
FEAT
Ours (no ref)
AFP
Targel Noise
gplearn-
x
0.0
FFX
0.001
DSR
0.01
ITEA
BSR
0.0
0.2
0.4
0.6
0.8
1.0
10a
104
10=
10a
10a
tradeoff's, even in presence of noise. Results are averaged over all 119 Feynman problems, for 10 random seeds and three target noises each as shown in the legend. The accuracy is computed as the fraction of problems for which the R² score on test examples is above 0.99. Models are ranked according to the accuracy averaged over all target noise.
Figure 5: Our model presents strong accuracy-ppee-ccopple
Conclusion
In this work, we introduced a competitive deep learning model for SR by using a novel numeric- symbolic approach. Through rigorous ablations, showed that predicting the constants in an expression not only improves performance compared to predicting a skeleton, but can also serve as an informed initial condition for solver to refine the value of the constants.
Inference single GPU for the forward pass of the Transformer.

---

Our model outperforms previous deep learning approaches by margin on SR benchmarks, and scales larger dimensions. Yet, the dimensions considered here remain moderate (D 10): adapting to the truly nigh-dimensional setup is an interesting future direction, and will likely require qualitative changes the data generation protocol. While our model narrows the gap between GP and DI. based SR, closing the gap also remains a challenge for future work.
This work opens up whole new range of applications for SR fields which require real-time inference. We hope that the methods presented here may also serve as a toolbox for many future applications of Transformers for symbolic tasks.
10

---

References
[1] Silviu-Marian Udrescu and Max Tegmark. Ai feynman: a physics-inspinee method for symbolic regression, 2020.
[2] M. Cranmer, Alvaro Sanchez- Sanchez-Gonzale Gonzalez. Peter W. Battaglia, Rui Xu, Kyle Cranmer, David N. Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. ArXiv, abs/2006. 1287, 2020.
[3] Jiří Kubalik, Erik Demerr, Jan Zegklitz, and Robert Babuška. Symbolic regression methods for reinforcement learning. IEEE Access, 9:139697-139711 2021.
[4] Erik Derner, Jirí Kubalík, Nicola Ancona, and Robert Babuška. Symbolic regression for constructing analytic models in reinforcement learning. ArXiv, abs/1903. 2019.
[5] Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. Symbolicgpt: A generative transformer model for symbolic regression. arXiv preprini arXiv:2106. 14131, 2021.
[6] Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Paras- candolo. Neural symbolic regression that scales, 2021.
[7] William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabricio Olivetti de Franca, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason H Moore. Contemporary symbolic regression methods and their relative performance. arXiv preprint (tti::2200.113351,, 2021.
[8] Nikos Aréchiga, Francine Chen, Yan- Ying Chen, Yanxia Zhang, Rumen Iliev, Heishiro Toyoda, and Kent Lyons. Accelerating understanding of scientific experiments with end end symbolic regression. ArXiv, abs/21 12,0402., 2021.
[9] Silviu-Marian Udrescu and Max Tegmark. Symbolic pregression: Discovering physical laws from raw distorted video. Physical review. E, 103 4- 4::43333,7, 2021.
[10] Anja Butter, Tilman Plehn, Nathalie Soybelman, and Johann Brehmer. Back to the formula lhe edition. 2021.
[11] Michael Schmidt and Hod Lipson. Age-fitness pareto optimization, In Genetic programming theory and practice VIII, pages 129-146. Springer, 2011.
[12] Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324(5923):81 2009.
[13] William La Cava, Tilak Raj Singh, James Taggart, Srinivas Suri, and Jason H Moore. Learn- ing concise representations for regression by evolving networks of trees. arXiv preprint arXiv:1807.0098) 2018.
[14] Trent McConaghy. Ffx: Programming Theory and Practice IX, pages 235-260. Springer, 2011.
Fast, scalable, deterministic symbolic regression technology. In Genetic
[15] 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-237, 2021.
[16] Fabricio Olivetti de França and Guilherme Seidyo Imai Aldeia. Interaction-traasffrmaatooo evolutionary algorithm for symbolic regression. Evolutionary computation, 29(3):367-390, 2021.
[17] Ignacio Arnaldo, Krzysztof Krawiec, and Una-May O'Reilly. Multiple regression genetic programming. In Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation, pages 879-886, 2014.
[18] Marco Virgolin, Tanja Alderliesten, and Peter A. Bosman. Linear scaling with and within semantic backpropagation--aasee based genetic programming for symbolic regression. In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO '19, page 1084-1092, New York. NY, USA, 2019. Association for Computing Machinery.

---

[19] Michael Kommenda, Bogdan Burlacu, Gabriel Kronberger, and Michael Affenzeller. Parameter identification for symbolic regression using nonlinear least squares. Genetic Programming and Eyolvable Machines, 21(3):471-501, 2020.
[20] Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne T Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXiv preprint arXiv: (9Xi::/912.04871,, 2019.
[21] Georg Martius and Christoph H Lampert. Extrapolation and learning equations. arXiv preprint arXiv:1610.02995 2016.
[22] Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. In International Conference on Machine Learning, pages 2018.
PMLR,
[23] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in neural information processing systems, pages 5998-6008, 2017.
[24] Guillaume Lample and François Charton. Deep learning for symbolic mathematics. arXiv preprint 1rXiv:191..014:1 2019.
[25] François Charton, Amaury Hayat, and Guillaume Lample. Learning advanced mathematical computations from examples. arXiv preprint arXiv:2006.0644 06462, 2020.
[26] François Charton. Linear algebra with transformers. arXiv preprint arXiv:2112.0088 2021.
[27] Christopher Hahn, Frederik Schmitt, Jens U Kreber, Markus N Rabe, and Bernd Finkbeiner. Teaching temporal logics to neural networks. arXiv preprint 20X1v22003.042/8, 2020.
[28] Stanislas Polu and Ilya Sutskever. Generative language modeling for automated theorem proving. arXiv preprint arXiv:2009.03:: 03393333,,,,, 2020.
[29] Stéphane d' Ascoli, Pierre-Alexanrre Kamienny, Guillaume Lample, and François Charton. Deep symbolic regression for recurrent sequences. arXiv preprint arXiv:2201.04660 2022.
[30] Nikos Arechiga, Francine Chen, Yan- Yan-Ying Chen, Yanxia Zhang, Rumen Iliev, Heishiro Toyoda, and Kent Lyons. Accelerating understanding of scientific experiments with end to end symbolic regression, 2021.
[31] Roger Guimerà, Ignasi Reichardt, Antoni Aguilar-Mogas, Francesco A Massueci, Manuel Miranda, Jordi Pallarès, and Marta Sales-Pardo. A bayesian machine scientist aid in the solution of challenging scientific problems. Science advances, (()))eav6971,
[32] Patrick Kidger. Sympytorch.
2021.
[33] Richard Zou Horace He. functorch: Jax-like composable function transforms for pytorch. https://gt thg.. cowp//yyooorrccc/ttttorrolo 2021.
[34] Sam Wiseman and Alexander M. Rush. Sequence-to-sequenec learning as beam-search opti- mization, 2016.
[35] Jerome H Friedman. Greedy function approximation: a gradient boosting machine. Annals of statistics, pages 1189-1232, 2001.
[36] Steven Strogatz. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering. Westview Press, 2000.
[37] Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression, 2020.
[38] T. Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P. Santiago, Daniel M. Fais- sol, and Brenden K. Petersen. Symbolic regression via neural-guided genetic programming population seeding, 2021.
12

---

[39] Sean Welleck, Peter West, Jize Cao, and Yejin Choi. Symbolic brittleness in sequence models: on systematic generalization symbolic mathematics. Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages
2022.
[40] Martin Vastl, Jonáš Kulhánek, Jirí Kubalík, Erik Derner, and Robert Babuška. Sym- former: End-to-end symbolic regression using transformer-basee architecture. arXiv preprint arXiv:2005.15764.222 2022. 2022.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A]
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either the supplemental material or as URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect the random seed after running ex- periments multiple times)? [No] We noticed the random seed that dictates model initialization had no influence on performance after few training epochs, therefore we trained our models on single random seed to avoid unnecessary computations, especially learning is costly.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes]
4. If you using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any assets either in the supplemental material or as URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you'l rourre using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable?? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
13


Reviews:
[{"review_id": "sWQ_Pe4tjFN", "reviewer": "Reviewer_LNce", "summary": "This paper proposed a transformer-based approach to perform end-to-end symbolic regression, that is, predicting simultaneously the function and the values of numerical constants. The trained model is augmented with several tricks at inference, including refinement, scaling, and bagging, etc. Experimental results show that the accuracy approaches to the state-of-the-art GP approach and the inference time is significantly reduced by several orders of magnitude.", "questions": "1. Could you make a further comparison with [29] instead of just mentioning “[29] infers one-dimensional recurrence relations”? It seems that both of them have the same architecture and similar tasks.\n\n2. Could you please provide more ablation studies to demonstrate the effects of: (1) mixture of distributions and different distribution shapes; (2) input points scaling at inference time.\n\n3. Could you please compare with other deep learning symbolic methods like [5] in experiments? It seems that [5] performs well on both running time and accuracy.", "limitations": "Please see questions listed above.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Symbolic regression is well studied with various machine learning technique. In contrast to addressing the problem in a two-step manner, the proposed approach directly predicts the function along with the value of constants. Such an end-to-end model is not surprisingly new. The key innovation is the employed tricks for improving the performance of trained model, and the employed techniques are technically sound.\n\nAs the accuracy is still lower than state-of-the-art GP approaches, the main benefit of the proposed approach is the reduced inference time. However, I am not sure whether inference time is the key factor of symbolic regression. It would be helpful to introduce more background on the importance of inference time.\n\nBesides, more comparison with other end-to-end symbolic regression approach and ablation studies would be quite helpful to demonstrate the effectiveness of the proposed approach.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "CMissa5zgmZ", "reviewer": "Reviewer_Hxn5", "summary": "This paper proposes a symbolic regression model using Transformers. It enables to pre-train the model using datasets from the other domain, which is helpful when the target domain dataset is relatively small. Moreover, this approach performs the parameter refinement with BFGS to further enhance the performance. Compared to the state-of-the-art approaches, which are mainly genetic programming based, this approach achieved not superior but competitive performance on benchmark tasks, and showed significantly smaller inference time.\n", "questions": "Please see the comment above.", "limitations": "Although the checklist says the code is available, I couldn't find the pointer to the code in the main text. It would always be nice to have a pointer to the code, even if it is included in the supplementary. \n\nAlgorithmic description of the proposed approach would help readers understand.", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Clarity\n\nThe motivation of each algorithmic component is described well. However, as I am not an expert in SR, it was hard for me to understand this paper in details. Probably it is because of my lack of the common knowledge in the SR community. For example, how is it pretrained? I couldn’t find any statement about the setting for the pretraining. Or is the method written in Section 2.1 for the pretraining? If so how is the post-training performed? The computational time written in 166 is for what? It should depends on the dataset, right? It was also not clear for me how the refinement of the parameter is performed in the end-to-end manner. \n\nSignificance\n\nBased on the authors’ statement, this approach is the first one that reaches competitive performance to GP-based approaches. Although the performance is not superior, the advantage of using pre-training over GP-based approach is promising.\n\nThe author’s claim the advantages of the proposed approach over the top performance approaches in its inference time. However, it was not evaluated whether this comes from the trick described in Section 2.2 or it is intrinsic to the approach using Transformer. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "1Bw8upTML-Z", "reviewer": "Reviewer_ye2k", "summary": "This paper provide a novel sequence prediction based solution to symbolic regression (SR). Compared to search-based SR methods, the proposed method is more efficient in inference and generalize to new tasks effortlessly. Compared to privious sequence based SR methods, the proposed E2E is able to scale to larger dimensions.", "questions": "One important issue for the sequence based SR is that people need to know to what extent the model is learned to overfit the pretraining data. The authors provided OOD evaluations in the exp section  and the dataset generation controlling params in the supples, which is helpful to this question. However, it would be better to 1) expliclitely show how different is the OOD symbolic samples differ from the pre-training data generation params? Better demonstrate this with concrete examples of OOD symbolic ground truth and  pretraining examples. And 2) it is helpful but not mandatory to: show that modifying the pre-training dataset distribution could positively contribute to but not greatly impact the OOD result.\n\n\nAs a new method targeting at sequence based SR, it is better to compare with prior works in this domain, which are symbolicGPT and NSRS. I am aware that official implementations seems to be not yet available, hence I'm not posting this as a must-have. However, since the E2E SR and symbolicGPT are similar in many ways, it would still be beneficial if we could see the simulated results: one way is to alter the E2E implementation so as to only predict skeletons, then apply coeff optimization to un-officially simulate the symbolicGPT.", "limitations": "This work does not explicitely show potential negative societal impact.\n", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nSymbolic regression is one challenging, theoretically and practically important task that incorperates pattern recognition, rule based reasoning, sequence modeling, discrete optimization, and beyond. \n\nCompared to GP based SR, the sequence prediction based SR scales to novel unseen problems effortlessly. Compared to previous works in the sequence prediction based SR such as symbolicGPT and NSRS, this work scales better, capable to deal with > 3 dimensional problems and is with better accuracy. This paper claims that using the proposed dataset generation scheme, the pre-trained E2E SR could generalize to both in-distribution and OOD problems.\n\nIn the formulation of sequence based SR, the most critical issue is how to sample the pre-training tasks such that the model learned on the generated samples could capture symbolic knowledge, while not overfit to the pre-training dataset. To address this, the authors provided a principled dataset generation workflow, and showed via experiments that E2E could generalize to OOD. The author also provided in the supples additional series of experiments on whether the E2E is memorize the data.\n\n\nThe authors provided implementations which helps the reproducibility.\n\nThough being a new SR paper that pushes the SOTA, this paper upscales the dimension from (<=3) only to (<10), as the author mentioned in line 287. Hence the method still remains improvable in scaling to larger dimensions.\n\nOther limitations dwell in the evaluation of the pre-training data generation procedure, which is critical to this sequence prediction type of SR formulation, but is under evaluated.\n\nI haven't checked the provided implementations carefully. If the results are fully reproducible and more evaluations could be given to pre-training data generation on OOD, based on my limited assessment, this paper could be given high merit.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept