ParFam -- (Neural Guided) Symbolic Regression via Continuous Global Optimization
Authors: Philipp Scholl
Katharina Bieker
Hillary Hauger
Gitta Kutyniok
Venue: ICLR 2025 Poster (2025)
Keywords: symbolic regression
continuous optimization
expressivity
transformers
supervised learning
Primary area: other topics in machine learning (i.e., none of the above)
TL;DR: We introduce the symbolic regression method ParFam, which achieves state-of-the-art by leveraging the structure of common physical laws, theoretically analyze its expressivity and extend it using a pre-trained transformer.

Abstract: The problem of symbolic regression (SR) arises in many different applications, such as identifying physical laws or deriving mathematical equations describing the behavior of financial markets from given data. Various methods exist to address the problem of SR, often based on genetic programming. However, these methods are usually complicated and involve various hyperparameters. In this paper, we present our new approach ParFam that utilizes parametric families of suitable symbolic functions to translate the discrete symbolic regression problem into a continuous one, resulting in a more straightforward setup compared to current state-of-the-art methods. In combination with a global optimizer, this approach results in a highly effective method to tackle the problem of SR. We theoretically analyze the expressivity of ParFam and demonstrate its performance with extensive numerical experiments based on the common SR benchmark suit SRBench, showing that we achieve state-of-the-art results. Moreover, we present an extension incorporating a pre-trained transformer network (DL-ParFam) to guide ParFam, accelerating the optimization process by up to two magnitudes. Our code and results can be found at https://github.com/Philipp238/parfam.

Published as a conference paper at ICLR 2025
GUIDED) SYMBOLIC REGRES- SION VIA CONTINUOUS GLOBAL OPTIMIZATION
Philipp Scholl, Katharina Bieker, Hillary Hauger & Gitta Kutyniok"
Mathematical Institute
LMU Munich
Munich, Germany
scholl@me mooll@mmahh 1mu de
ABSTRACT
The problem of symbolic regression (SR) arises in many different applications, such as identifying physical laws or deriving mathematical equations describ- ing the behavior of financial markets from given data. In this paper, present our new approach ParFam that utilizes parametric families of suitable symbolic functions translate the discrete symbolic regression problem into continuous one, resulting in a more straightforward setup compared to current state-of-the- art methods. In combination with global optimizer, this approach results in a highly effective method to tackle the problem of SR. We theoretically analyze the expressivity of ParFam and demonstrate performance with extensive numerical experiments based on the common SR benchmark suit SRBench, showing that we achieve state-of-the results. Moreover, we present an extension incorporating a pre-trained transformer network (DI-PPaFaam) to guide ParFam, accelerating the optimization process by up to two magnitudes. Our code and results can be found at https:/ggttabbccom//hiiig///////
INTRODUCTION
Symbolic regression (SR) aims to discover concise and interpretable mathematical functions that accurately model input-output relationships. This focus on simplicity crucial for applications requiring model analysis and trustworthiness, such in physical chemical sciences (Quade al., 2016 Angelis et al., [023) Wang et al. 2019). SR finds broad application in diverse fields, including ecosystem dynamics (Chen al. 2019), solar power forecasting (Quade et al. 2016). financial
market analysis (Liu & Guo 2023) , materials science (He & Zhang 2021). and robotics (Oplatkoval
& Zelinka 2007). The growing body of SR research, as evidenced by the increasing number of publications (Angelis et al., (203)), underscores its significance.
SR is a regression task in machine learning that aims to find an accurate model without any assump- tions by the user related to the specific data set. Formally, a symbolic function f R + R that accurately fits a given data set (I,, 1...=.....N S IR" for all data points, in the case of noise, yi f(x;) for all E (1,..., N). Unlike other regression tasks, SR aims finding a simple symbolic and thus interpretable formula while assuming as little as possible about the unknown function. In contrast to SR, solutions derived via neural networks (NNs), for instance, lack interpretability. Traditional regression tasks, on the other hand, typically assume strong structure of the unknown function, such as lincarity or polynomial.
sought, i.e., it should satisfy yı = (a;)
To tackle SR problems, the most established methods are based on genetic programming (Augusto & Barbosa 2000 Schmidt & Lipson| 2009 2010 Cranmer 2023) and nowadays many algorithms
incorporate neural networks (Martius & Lampert 2017) Udrescu & Tegmark 2020; [esai & Stra-
chan 2021 Makke al.) (2022). However, despite significant effort, many methods struggle to
consistently find accurate solutions on challenging benchmarks. La Cava et al. (2021) evaluate 13 SR algorithms on the SRBench ground-truth problems: the Feynman (Udrescu Tegmark, 2020)
*P. S. and G. K. are also affiliated with the Munich Center for Machine Learning (MCML), Germany. G. K. is further affiliated with the German Aerospace Center (DLR), Institute of Robotics and Mechatronics, and University of Tromso, Norway.

---

Published as a conference paper at ICLR 2025
and Strogatz (La Cava et al. 2016) problem sets. Both data sets consist of physical formulas with varying complexities, where the first one encompasses 115 formulas and the latter 14 ordinary dif- ferential equations. Most algorithms achieved success rates below 30% on both datasets within an 8-hour time limit, with only Al Feynman (Udrescu & TTegmark 2020) showing better performance. Furthermore, these results deteriorate significantly in the presence of noise LLa Cava et al. 202T Cranmer )023
In this paper, we introduce ParFam, novel SR algorithm that leverages the inherent structure of physical formulas. By translating the discrete SR problem into a continuous optimization problem, ParFam enables precise control over the search space and facilitates the use of gradient- opti- mization techniques like basin-hopping (Wales & Doye 1997). While ParFam is not the first method to employ continuous optimization for symbolic regression, it aims to enhance the translation the continuous space. By doing so, ParFam becomes the first SR method based on continuous opti- mization to achieve state-of-thh-ar performance. To accelerate ParFam at the cost of flexibility, we extend it to DL-ParFam which incorporates pre-trained Set Transformer (Lee et al. 22119) to guide the optimization process for ParFam. ParFam is detailed Section 2.T its expressivity is analyzed in Section 2.2] and DL-ParFam is introduced in Section 2.3 Notably, despite its simplicity, ParFam achieves state-of-the of-t-----tttttttt the results on the Feynman and Strogatz datasets, as demonstrated in Section
Our Contributions Our key contributions are as follows:
1. Introduction of ParFam, novel method for SR improving compared to existing continuous optimization-buss based SR algorithm by leveraging the inherent structure of physical formulas and the expressivity of rational functions to translate SR into efficiently solvable con- tinuous optimization problem, by avoiding the need for nested basis functions. This results in the following advantages: (1) Enabling gradient-based optimization techniques while avoiding exploding gradients, (2) enhanced interpretability, and (3) efficient but simple and user-friendly setup.
2. Thorough theoretical analysis of the expressivity of ParFam, showing its high expressivity despite its pre-defined structure necessary for continuous optimization.
3. Introduction of DI-ParFam, an extension of ParFam based pre-trained Set Transformer (Lee et al. 2019),, which guides ParFam and, therefore, accelerates its search by up to 100 times.
4. Extensive benchmarks showing state-of-the-an performance of ParFam and DL-ParFam and significantly better results than other methods based on continuous optimization.
Related work Most SR algorithms approach the problem in two steps. First, they search for the analytic form of the target function in the discrete space of functions and then optimize the coefficients via continuous optimization techniques like BFGS (Nocedal & Wright 2006).
Traditionally, genetic programming was used to heuristically search the space of equations given some base functions and operations (Augusto & Barbosal 2000; (Shhmitt & Lipson 2009: 2010
Cranmer 2023). However, due to the accomplishments of NNs across diverse domains, numerous researchers aimed leverage their capabilities within the realm of SR. [ddrescu & Tegmark 22220) for instance, have employed an auxiliary NN to evaluate data characteristics.
Petersen al. (222T) rely on reinforcement learning (RL) explore the function space, where a policy, modeled by a recurrent neural network, generates candidate solutions. Mundhenk et al. (202T) combined this concept with genetic programming such that the RL algorithm iteratively learns to identify a good initial population for the GP algorithm, resulting in superior performance compared individual RI. and GP approaches. Similarly, (Sun al.) (2022) rely on Monte Carlo tree search to search the space of expression trees for the correct equations.
Inspired by the success of pre-training of models on large data sets other machine learning tasks KKallnn et al. 2020 Devlin et al. 20188 Brown et al.) 2020; Chen et al. 2020) and for mathematical data (Lample & Charton| 2019), there have been many attempts in recent years to leverage pre-
training for SR. Biggio al. ( I build upon the data generation from Lample & Charton (219) to train a neural network to predict the skeleton of a function symbolically, the constants of the skeleton
then found via BFGS. Kamienny et al. (022) extend this approach predict the skeletons
2

---

Published as a conference paper at ICLR 2025
directly and only use BFGS for fine-tuning afterwards. training with a combination of RL (Petersen et al., 2021). GP (Mundhenk et al. 2021). Al Feynman (Udrescu & Tegmark 2020) and linear models.
Landajuela al. (2022) incorporate pre-
In contrast, only a few algorithms, such as FFX McConaghy 2011) and SINDY (Brunton et al, 2016), share ParFam's 's approach of merging the search for the analytical form with coefficient op- timization. These methods utilize a model with linear parameters, enabling efficient coefficient estimation via sparse linear regression. To expand the search space, they generate large set of features by applying base functions to the input variables. However, this linear parameterization restricts the search space, as cannot model non-linear parameters within the base functions.
The closest method to ParFam is EQL with division (Martius & Lampert 2017 Sahoo et al. 2018). which makes use of continuous optimization and overcomes the limitation of FFX and SINDy by utilizing small NNs with sin, cos, and the multiplication as activation functions. The goal of EQL to find sparse weights such that the NN reduces an interpretable formula. However, while EQL applies linear layers between the base functions, ParFam applies rational layers. Thereby, EQL usually needs multiple layers to represent the most relevant functions, which introduces many redundancies, inflates the number of parameters, and complicates the optimization process. More- over, EQL relies on the local minimizer ADAM (Kingma & Ba 2014) for coefficient optimization. On the contrary, ParFam leverages the reduced dimensionality of the parameter space by applying global optimization techniques for the parameter search, which mitigates the issues of local minima. Furthermore, ParFam maintains versatility, allowing for the straightforward inclusion of different base functions, while EQL cannot handle, c.g., the exponential, logarithm, root, and division within unary operators. In recent years, several extensions of EQL and similar approaches have been pro- posed. DySymNet (Li et al., 2024) employs a structure akin to EQL but optimizes the architecture through reinforcement learning. MetaSymNNt (Li al., 2023) builds on the EQL framework by incorporating evolvable activation functions and rules for dynamically modifying the architecture during training. Similarly, Dong et al. (2024) enhance the classical EQL. network by introducing an activation function that is evolved using genetic programming.
2 METHODS
In the following section, we first introduce ParFam novel approach, that exploits well-suited representation of possible symbolic functions to which an efficient global optimizer can be applied. Afterwards, theoretically analyze the expressivity of the considered parametric families, showing that ParFam is quite expressive despite the restrictions due to the predefined structure. We then introduce DL-ParFam as an extension of ParFam which incorporates pre-training on synthetic data to inform the parameter choices of ParFam and, therefore, speed up the learning process.
2.1 PARFAM
The aim of SR is to find a simple and thus interpretable function that describes the mapping under- lying the data (xi, y1) 1..... without additional assumptions. Typically, set of base functions, such as (+,-,- exp, sin, /.. is predetermined. The primary goal of an SR algorithm is to find the simplest function that uses only these base functions to represent the data, where simplicity is usually defined as the number of operations. In Subsection 2.1.1 we introduce the architecture of ParFam and in Subsection 2.1.2 we explain its optimization.
2.1.1 THE STRUCTURE OF THE PARAMETRIC FAMILY
The main goal of ParFam is to translate SR into continuous optimization problem to enable the use of gradient and higher-order derivative information to accelerate the search. achieve this, we construct parametric family of functions in a neural network-liks structure presented in Figure [ where ...... Qx+1 are rational functions and J1..... Ss are the unary base functions, which cannot be expressed rational functions, like sin, v exp, etc. The difference from standard residual neural network with one hidden layer is that we use rational functions instead of linear connections between the layers. Furthermore, we apply physically relevant activation functions 9I....9k which may differ in each neuron.
3

---

Published as a conference paper at ICLR 2025
311
91
Q
Qu1
22
M
() (
Q:(3
Figure The architecture of ParFam: ParFam can interpreted as a residual neural network with one hidden layer. Instead of linear weights between the layers, it applies rational functions Qi(+) = Pa (.)...)). Furthermore, the standard basis functions are substituted by physically relevant functions like sin, exp, v. etc. The learnable parameters are the coefficients of pay and Pat-
Since the network has only one hidden layer, we can write it in a compact form as (o(z)=Qx+1l (z, 9i (Q1(z)),
99:..9gx( (x(()(z))),
(1)
where z E IR the input vector. The learnable parameters 0 R are the coefficients of the polynomials, i.e., of the numerators and denominators of Q.....,Qu Qk+1- The degrees df and df. (1....,k+)) k+1], of the numerator and denominator polynomials of QI.... QK+1, respectively, and the base functions 91...., Is are chosen by the user. Depending on the application, custom functions can be added to the set base functions. This versatility and its simplicity make ParFam highly user-friendly tool, adaptable to a wide range of problem domains. In Appendix A we explain how to incorporate specific base functions to avoid numerical issues and further implementation details.
It is possible to extend the architccture shown in Figure to multiple layers to arbitrary func- tions. However, the main motivation for the proposed architecture that it consists of a single hidden layer, due the high approximation qualities of rational functions and the general structure of common physical laws (Woan 2000). Employing a single hidden layer offers several advantages: it reduces the number of parameters, simplifies optimization by mitigating issues such explod- ing or vanishing gradients caused by nested functions, and enhances interpretability since it avoids complicated and uncommon compositions such as sin cos (Woan) 2000), which many algorithms enforce to avoid as well (Petersen et al. 2021 Landajue et al. 2022). In Section [.2 we analyze the expressivity of our architecture (using one hidden layer) and in Section 3 we show that the structure is not only flexible enough to recover many formulas exactly, but also has the best approximation capabilities among all tested algorithms.
2.1.2 OPTIMIZATION
The goal of the optimization is to find the coefficients of the rational functions Q1,..., 1 such that fe approximates the given data (z,, Ii..= thus, we aim to minimize the mean squared error (MSE) between y/ and fo;).. As we aim for preferably simple functions derive interpretable and easy-to-ane tonallyee results, regularization term R(0) is added to encourage sparse parameters. In total, we consider the loss function
fe(2)))2 (( = E -olc))) AR(O).
where > 0 is a hyperparameter to control the weight of the regularization. Here, we choose R(O) |0|| as surrogate for the number of non-zero parameters, which is known enforce sparsity Bishhop 2006, Goodfellow et al. 2016). In Appendix A we discuss how to deal with the regularization of the coefficients of rational functions detail, to ensure that the regularization is applied effectively and cannot be circumvented during optimization.
Although the SR problem now transformed into continuous optimization problem, due to the presence of many local minima, it is not sufficient to apply purely local optimization algorithms like gradient descent or BFGS (Nocedal Wright 2006). This is also shown in our comparison study in Appendix B To overcome these local minima, we instead rely on established (stochastic) global optimization methods. Here, we choose the so-called basin-hopping algorithm, originally introduced by Wales & Doye (1997). which combines local minimizer, e.g., BFGS (Nocedal & Wright| 2006). with a global search technique inspired by Monte- Carlo minimization (Li & Scher- aga 1987) to cover larger part of the parameter space. More precisely, the implementation
4

---

Published as a conference paper at ICLR 2025
provided by the SciPy library (Virtanen al.) 2020). The basic idea of the algorithm is to divide the complex landscape of the loss function into multiple areas, leading to different optima. These the so-called basins. The random perturbation of the parameters allows for hopping between these basins and the local search (based on the real loss function) inbetween improves the results and ensures that global minimum is reached if the correct basin is chosen. For the acceptance test, the criterion introduced by Metropolis et al. (1953) is taken.
Following the optimization with basin-hopping, finetuning routine is initiated. In this process, co- efficients that fall below a certain threshold are set to O, and the remaining coefficients are optimized using the L-BFGS method, starting from the previously found parameters. The threshold grad- ually increased from 10-$ to 10 to encourage further sparsity in the discovered solutions. This step has been found to crucial in enhancing the parameters initially found by basin- masin-hopping.
2.2 EXPRESSIVITY OF PARFAM
Since the structure of the parametric family in ParFam restricts the search space, it is interesting to investigate the expressivity of our approach in detail. For this, we aim to quantify the ratio of the number of functions of complexity / N that can be represented using ParFam the number of functions with the same complexity that cannot be represented. We follow the modeling approach by [ample& & Charton 2 who used expression trees represent functions and investigate the number of expression trees, which is a common structure used in genetic algorithms. Examples of expression trees are given in Figure 2 Note that this approach does not take into account that some trees model the same expression, e.g., 21 + 22 and 22 + 21 will refer to different trees.
C00
DOS
21
#2
21
(b) Unary-binary expres- sion tree with inter- nal nodes, representable by ParFam, for exp(21) + cos(22).
(a) Binary expression tree with 2 internal nodes for 4:2
(c) Unary-binary expres- sion tree with 2 internal nodes, not representable by ParFam, for cos( /222).
Figure 2: Examples for the different kinds of expression trees counted by bt, ct. and dj.
We call node without children leaf, node with child unary node, and node with two children a binary node. The nodes with children, i.e., the unary and binary nodes, are referred to as internal nodes. We call a tree binary tree if all its internal nodes are binary, and unary-binary tree if all its internal nodes are either unary or binary. We assume that there are n different possible leaves (the variables 21..... 2,), & different unary nodes (the non-rational functions 91...., Ik),, and different binary nodes (the binary operations +, -/..;; can be set to 4 in general). We do not model constants separately, as these can be incorporated either as part of each node as additional leaves.
To quantify the expressivity of ParFam, we define as the number of unary-binary trees with / internal nodes that can be represented by ParFam and dy as the number of all unary-binary trees with internal nodes (including the trees that cannot be represented by ParFam). Examples for these different types trees are shown in Figure 2 The tree Figure 2c cannot be modeled by ParFam since there is path from the root (cos) to leaf (I22,, that contains more than unary node. Such paths represent compositions explained in Section 2. Typical formulas from the Feynman dataset Udrescu & Tegmark (020 for instance, have a complexity of around 10. For example, the formula m sis sin(n0 (n#/ sin (Feynman 1.30.3) has complexity of
unary functions 91, 9k which are omitted by ParFam as
Our goal in this section is to compute an estimate for the ratio c//h. The proofs rely mostly on the idea of generating functions (Will 2005) and can be found in Appendix B together with additional context. We start by stating an approximation to proven in Appendix E.5
5

---

Published as a conference paper at ICLR 2025
Theorem 2.1. For (ci)LeN, the number of unary-binary trees expressible by ParFam with complexity I. it holds that
G
2bz
(vo VAR(++1
8V4m(1=1)5
-1
UI 4V/(U+1/5 +0
(3)
with some constants Vo, U1 R depending on the number of binary operators b, number of unary operators k, and number of variables
In Appendix E.5 we additionally compute the exact formulas for to and t1- Moreover, we show the approximations and the true values cj in Figure Sal revealing that already the first-order ap- proximation of e is quite close to the exact one. Since we are interested in a/d,, we also need an approximation di
Theorem 2.2. For (di)ters the number of unary-binary trees with complexity I, holds that
d1
(V1-
Vaml++aps 8V/x((+1)))
3
V +06TT(1-712 (4)
Srg
2br 4T
8V1-
where 11,2 k+22bnd:22 bbbbkbbbbb/ and k and n denote the number of binary operators, the number of unary operators, and the number of variables, respectively.
122
The proof of Theorem 2.2 is given in Appendix E.7 Theorem and 2.2 yield an approximation for the expressivity of ParFam q/ dr. To simplify the approximations, we disregard the constant and polynomial terms in / yielding //d; z (f2/x1)'. Table shows r2/x1 for b = 4 and varying values for lt and n. The ratio is mainly between 0.9 and 1.0, revealing that especially for formulas of low complexity, i.e., small 1, is relatively close to di, showing the high expressivity of ParFam despite its restrictions. For example, for n = and 3, Table yields f2 = 0.9799. Therefore, d; 2 0.9799 holds. E.g, for I = 5, ParFam covers 90.25% of formulas, and, for = 10, 81.62%.
Table I: Approximation of G41/a ~
given by Theorem 2.1 and 2.2 for = 4 and varying
411/4
values for k and
n k:
2
3
4
5
6
0.9712
0.9356
0.9020
0.8713
0.8435
0.8183
0,9881
0,9712
0.9533
0.9356
0.9185
0.9020
0.9931
0.9827
0.9712
0.9593
0,9474
0,9356
0,9954
0.9881
0.9799
0.9712
0.9623
0.9533
0,9966
0,9912
0.9849
0.9782
0.9712
0.9641
0.9974
0.9931
0.9881
0.9827
0.9770
0,9712
0.9979
0.9944
0.9903
0.9859
0.9811
0,9762
0.9983
0.9954
0.9919
0.9881
0,9841
0,9799
0,9985
0.9961
0.9931
0.9899
0.9864
0.9827
2.3 DL-PARFAM
The choice of model parameters for fa, specifically the basis functions I1...., and the degrees of the polynomials Q QI...., Qk+1, is crucial yet challenging. A highly general parametric family can lead to a complex optimization problem, while restrictive family may hinder the discovery of the correct expression. If the time constraint allows it, iterating through different model parameters can be an effective solution for this problem, as demonstrated in our experiments in Section
To solve this problem with tighter time constraints, we introduce DI-ParFam: A neural-guided version of ParFam, which follows the recent emergence of pre-trained neural networks for symbolic regression (Biggio et al.) 2020 202 Kamienny et al. 2022 Holt et al., (023) to predict the correct model parameters for the given data. The idea behind these approaches is to leverage the simple generation of synthetic data for symbolic regression, to train networks to map data (I;, ys)i=1.....N directly- or indirectly using further optimization schemes- to the corresponding solution. After a costly pre-training step that can be done offline and only has to be performed once, this acquired knowledge be used speed up the learning process considerably.

---

Published as a conference paper at ICLR 2025
DL-ParFam
Model Parameters
Degrass numeralor
funclions
d dlpp
(It L.
girl
ParFam
6"
fo : -
Degrees denominator
Figure 3: DL-ParFam first applies the pre-trained neural network input data (æ;; :1..=...... which outputs the model parameters for ParFam: the degrees of the polynomials used in Q1,.., Qk+1 and the basis functions 91, Sk- Afterwards, ParFam can run using these settings to find the best parameters 0y and, therefore, identify the best fitting function fa.
However, pre-training is a highly complicated task, because of the complex data distribution and the need to be able to handle flexible data sets in high dimensions (Biggio et al. 2021 Kamienny et al.) 2022). Furthermore, networks trained on the symbolic representation of functions fail to incorporate invariances in the function space during training, e.g., x+ and y I are seen as different functions, pointed out by Holt al. (2023), which possibly complicates the training. Holt al. (2023) solve this by evaluating the generated function to compute the loss and update the network using RL. However, evaluating each function during the training instead of comparing its symbolic expression with the ground truth is computationally expensive. Moreover, due to the non-differentiability the network has be optimized using suitable algorithms like policy gradient methods. The approach most similar to DL-ParFam is SNR (Liu et al. 2023), which uses a pre-trained SET-Transformer to predict a mask for active connections within SymNet-a symbolic neural network similar to EQL. During inference, these predictions are further fine-tuned using RL. With DL-PPrFam, see Figure we aim to combine the best of both worlds, by predicting the model parameters for ParFam which are more robust to functional invariances than symbolic representa- tions. This approach yields significant improvements compared to ParFam and other pre-training- based SR methods:
DL-ParFam strongly reduces the computational burden of ParFam by predicting model parameters.
DL-ParFam predicts the structure of the function, which can be directly compared with the ground truth and, thereby, avoids the evaluation of the predicted function on the data in every training step and yields an end-to-end differentiable pipeline, while being able to manage invariances in the function space.
The idea behind DL-ParFam that we train a neural network on synthetic functions to learn to predict the correct model parameters for ParFam from the data (I;, U:)i=1.....- After the training is finished, the network is used on real data to predict the model parameters, with which ParFam then aims compute the underlying function. We implement the neural network for DL-ParFam using a Set Transformer (Lee et al. 2019) as the encoder and fully connected ReLU network as the classifier. To train it we generate synthetic functions along with the corresponding data points (I;, 3s)i=1....N and one-hot encoded model parameters. Note that we do not follow the data gener- ation introduced by Lample Charton| (2019) as other approaches (Kamienny al., 2022 Biggio et al., 2021 Holt et 2003) do since we require direct access to the model parameters for ParFam, which would be complicated to extract from the formulas generated by Lample Charton 2019).
The x; are sampled from uniform distribution
[1, 5]", where 1 < < 9. The functions
along with f(x,) y:)i=1.....N are generated as described in Appendix C The input data to the
neural network is normalized to be between -1 and improve the training and make the model applicable to different input ranges. The network is trained by optimizing the cross-entropy between its output and the model parameters using the Adam optimizer with learning rate of 0.0001 and
7

---

Published as a conference paper at ICLR 2025
step learning rate scheduler. Further implementation details and the hyper- perer-errmmeers necessary to replicate the experiments are provided in Appendix D
The network was trained about 4,000,000 functions for 88 epochs after it stopped improving for epochs. The creation of the data took 1.3h and the training took 93h on one RTX3090. We further fine-tuned the network noisy data on 4,000,000 noisy functions for 30h one RTX3090. For comparison, end-to-end (E2E) (Kamienny et al.) (002)) made use of 800 GPUh and NeSymRes (Big-) gio al., 2021), which is only trained for problems up dimension 3, used 10,000,000 functions and was trained on a GeForce RTX 2080 GPU for 72h.
3 BENCHMARK
After the introduction of the SR benchmark (SRBench) by La Cava et al. (2021), several researchers have reported their findings on SRBench's 's ground-truth and black box data sets, due to the usage of real-world equations and data, the variety of formulas, the size of the benchmark, and comparability to competitors. The ground-truth datasets are synthetic datasets following real physical formulas and the black-box datasets real-world datasets for which the underlying formulas are unknown. These data sets are described more detail in Appendix F In this section, we evaluate ParFam and DL-ParFam on the SRBench data sets and report their performance in terms of the symbolic solution rate, the coefficient of determination R², and their training time demonstrating the strong performance of ParFam and the additional speed up achieved by DL-ParFam.
SRBench ground-truth problems
Symbolic Solution Hata
Accuracy Solution nate
Training time (5)
Pysk
Terget Numiee
ParFam
X
XB
0,0
0.01
LCG2
0.1
APP.FE
ooleam
ITEA
0peion
EndtuBrd
Figure 4: Mcan results on the SRBench ground-truth problems. Following SRBench terminology, training time refers the time each algorithm requires compute result for specific problem, which corresponds to inference time for pre-trained methods.
Competitors We include the results reported by SRBench for 14 SR algorithms and extend them by running EndToEEd (Kamienny al. 2022), uDSR (Landajuela al. 2022). and PySR (Cranmer) 2023) on our machines. Further information on the algorithms and the chosen hyper- parameters can be found in Appendix
hyper-parameters
To ensure comparability with the results evaluated on SRBench, ation metrics as La Cava et al. (2021). For the ground-trute problems, we first report the symbolic solution rate, which is the percentage of equations recovered by an algorithm. Second, we consider the coefficient of determination
Metrics
the same evalu-
i= -u -
R =1
(5)
ΣN (1 7)
where Ú; = (I;) represents the model's prediction and the mean of the output data y. The closer R² is 1, the better the model describes the variation in the data. It is widely used measure for
8

---

Published as a conference paper at ICLR 2025
R²
Treining
12y3R
*Operon
-SBP-GP
GP GOMEA
XGB
TEPLEX
Rondomforest
TTEA
Kernel Ridge
AdaBoost
*AFP
"AFP
"gpleamn
MLP
+FXX
1D5R
*dDSR
*MRGP
+BSR
*AFeymman
-0.25 0.00 0.25
0.50
0.75
1.00
10"
10a
10a
10r
10P
10"
20
10"
Figure 5: Median H², formula complexity, and training time on the 77 black-box problems from SRBench (La Cava et al., 2021) with at most 10 independent variables. The asterisk indicates that it a symbolic regression method.
goodness-of-He since it is independent of the scale and variation of the data. Following La Cava et al. (021) we report the accuracy solution rate for each algorithm, defined as the percentage of functions such that R > 0.999. The original data sets do not include any noise. However, similar to La Cava al. (2021), we additionally perform experiments with noise by adding E N N(0, o² N Ei= y to the targets Yi, where denotes the noise level. For the black-box problems we report the median R² and the median complexity of the formula, as defined in La Cava et al. (2021). b, since the symbolic solution rate is not defined in this case.
Hyperparameters The hyperparameters of ParFam can be divided into two subsets. The first subset defines the parametric family (fe)oeRm, e.g., the degree of the polynomials and the set of base functions. A good choice for this set is highly problem-dependent However, in the absence of prior knowledge, it is advantageous to select a parametric family that is sufficiently expansive to encompass wide range of potential functions. In this context, we opt for sin, exp, and / as our base functions. For the first layer rationals QI...., Qu. we set the degrees of the numerator and denominator polynomials to 2. For Qk+1, we set the degree of the numerator polynomial to 4 and the denominator polynomial to 3. This choice results a parametric family with hundreds of pa- rameters, making challenging for global optimization. To address this issue, iterate for ParFam through various smaller parametric families, each contained in this larger family, see Appendix H for details. The second set of hyperparameters defines the optimization scheme. Here, we the regularization parameter to À = 0.001, the number of iterations for basin-hopping to 10, and the maximal number of BFGS steps for the local search 100 times the dimension of the problem. Our choice of parameters for ParFam and DL-ParFam are summarized in Table 7 in Appendix I The hyperparameters for the pre-training of DL-ParFam can be found in Appendix D
Results Following La Cava et al. (2021) we allow a maximal training time of 8 CPU hours and maximal number of function evaluations 1, 000, 000 on the ground-truth data sets and 48 CPU hours on the black-box problems. In Figure we present the mean of the symbolic solution rate, the accuracy solution, and the training time on both data sets together. PySR, ParFam, AI Feynman, DL-ParFam, and uDSR outperform all other competitors by substantial margin (over 20%) when it comes to symbolic solution rate. Among those algorithms, PySR performs the best, followed by ParFam. These two algorithms are also the most robust noise, where it is important to notice that PySR is the only method that incorporates a noise filter (Cranmer) 2023).
Concerning the accuracy solution, ParFam outperforms all competitors with and without noise, fol- lowed by PySR, MRGP, Operon, uDSR, and DL-ParFam. These two metrics underscore the strongly competitive performance of ParFam with the current state-of-the-art While DL-ParFam performed slightly worse, beats most of the established methods in both metrics, while being up to a hun- dred times faster than its competitors, which was the goal of incorporating pre-training into ParFam.

---

Published as a conference paper at ICLR 2025
However, DL-ParFam's 's ability to recover the symbolic solution is notably hindered under low-noise conditions. The only other pre-training-bar baased method, EndTooEnd, performed worse all three met- rics. We performed the experiments without tuning the hyperparameter X. To assess the sensitivity of the results with respect to X, see Appendix K
In Figure 5 we present the median R², formula complexity, and training time across the 77 black- box problems with a maximum of 10 independent variables. While ParFam's performance is slightly weaker compared to the physics datasets discussed earlier, it remains among the rop-performing algorithms for these real-world datasets. Results of DL-ParFFa on the black-box data sets with at most variables (the limit of the current DL-ParFam) are shown in Appendix [
Due to the proximity of EQL (Martius & Lampert 2017; Sahoo et al., 2018) and ParFam, we deem a comparison between these two methods as highly interesting. However, the restricted expressivity of EQL makes it an unfair comparison on the whole Feynman and Strogatz dataset. For this reason, we show the results for EQL on reduced benchmark in Appendix|M] MM] To compare DL-ParFam with fur- ther pre-training-baser methods, we show the results of DL-ParFam and NeSymRes (Biggio et al.) 2021] on a reduced subset of Feynman and Strogatz in Appendix NN since NeSymRes cannot han- ale expressions with more than 3 variables. We further compare DL-ParFam with model-parament selection guided by Bayesian optimization in Appendix [ For results ParFam on the Nguyen benchmark {Uy et al. 201T) and comparisons with algorithms that were not tested SRBench, like SPI. (Sun et al. 2022) and NGGP (Mundhenk et al., (021). see Appendix P To show the robustness of ParFam and, especially, DL-ParFam to different data domains we refer to the SRSDD-Fyymman benchmarks (Matsubara et al. 2024) in Appendix Q
4 DISCUSSION AND CONCLUSION
This work introduces ParFam, accompanied by a theoretical analysis and extensive experiments. ParFam is the first continuous optimization-bartt based SR algorithm to match the performance of top genetic programming methods. Furthermore, present a novel pre-training approach that sig- nificantly outperforms existing methods and offers a substantial speed advantage over traditional competitors.
Limitations The parametric structure of ParFam is its greatest asset in tackling SR but also its main constraint, as it limits the function space. However, in Section 2.2 we theoretically prove that this limitation not severe. Section further demonstrates that algorithms theoretically capable of identifying specific formulas often fail in practice, while ParFam, despite its constraints, still finds highly accurate approximations. Another drawback is the computational expense of solving high- dimensional problems (>10 variables) with global optimizer, as the number of parameters grows in O where n is the number of variables and the polynomial degree. For DL-ParFam, the biggest challenge is the costly pre-training of the transformer network, making it less flexible than ParFam. However, since training is done offline with synthetic data, it can be reused for various SR tasks and its pre-training is faster than other pre-training- based methods.
Future Work Several avenues of ParFam and DL-ParFam remain unexplored, encompassing diverse forms of regularization, alternative parametrizations, and the potential incorporation of custom-tailored optimization techniques.
Acknowledgement This work of G. K. was supported in part by the Konrad Zuse School of Ex- cellence in Reliable AI (DAAD), the Munich Center for Machine Learning (BMBF) as well as the German Research Foundation under Grants PPG-98 KU 1446/31-1 and KU 1446/32-1. Fur- thermore, P.S. and G. K. acknowledge support from LMUxxcellent, funded by the Federal Ministry of Education and Research (BMBF) and the Free State of Bavaria under the Excellence Strategy of the Federal Government and the Länder well by the Hightech Agenda Bavaria and additional support by the project "Genius Robot" (011S24083), funded by the Federal Ministry of Education and Research (BMBF). Moreover, this work of H. H. is supported by the DAAD programme Konrad Zuse Schools of Excellence in Artificial Intelligence, sponsored by the Federal Ministry of Educa- tion and Research.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Dimitrios Angelis, Filippos Sofos, and Theodoros E. Karakasidis, Artificial intelligence in physical sciences: Symbolic regression trends and perspectives. Archives of Computational Methods in Engineering, (0(6))3845-386) 2023. doi: 101007/11831-003--000
D. Adriano Augusto and Helio J.C. Barbosa. Symbolic regression via genetic programming. In Proceedings. Vol.1. Sixth Brazilian Symposium on Neural Networks, pp. 173-178, 2000. doi: (0000///BB.22000088880
Deaglan J. Bartlett, Harry Desmond, and Pedro G. Ferreira. Exhaustive symbolic regression. IEEE Transactions Evolutionary Computation, pp- 1-1, 2023. doi: 10. 109/1EVC..22333220055
Luca Biggio, Tommaso Bendinelli, Aurelien Lucchi, and Giambattista Parascandolo. A seq2seq approach to symbolic regression. In Learning Meets Combinatorial Algorithms at NeurIPS2020, 2020.
I Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurélien Lucchi, and Giambattista Parascan- dolo. Neural symbolic regression that scales. Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learning Research, pp. 936-945, PMLR, 2021. URL http:/h/pppoceeednngeeoooooo//ee mmrr..mrr...... wttpppppppp.////////l/g/g....... att.mommlla.mmllllll.
Christopher M. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer New York, 2006.
Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in neural information processing systems, 33:1877-1901, 2020.
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, 113(155:39322-337, 2016. doi: 10. 1..77//paa..5113841111
Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. A simple framework for contrastive learning visual representations. International conference on machine learning, pp. 1597-1607. PMLR, 2020.
Yize Chen, Marco Tulio Angulo, and Yang- Yu Liu. Revealing complex ecological dynamics via symbolic regression. BioEssays, 41 (12)::190006) 2019. doi: 10. 00.100//bess.211900066
Miles Cranmer. Interpretable machine learning for science with pysr and symbolicregression. jl. arXiv preprint arXiv:2305.0/558 2023.
Saaketh Desai and Alejandro Strachan. Parsimonious neural networks learn interpretable physical laws. Scientific reports, 1005s111110 111):12761,
Jacob Devlin, Ming- Ming Chang. Kenton Lee, and Kristina Toutanova. Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv: 1810.04805 2018.
Junlan Dong, Jinghui Zhong, Wei-Li Liu, and Jun Zhang. Evolving equation learner for symbolic regression. IEEE Transactions on Evolutionary Computation. 2024.
Ian J. Goodfellow, Yoshua Bengio, and Aaron C. Courville. Deep Learning. Adaptive computation and machine learning. MIT Press, 2016. URL http://wwwtl wwwpweeweeepeeeiingpoooo...o org
Mu He and Lei Zhang. Machine learning and symbolic regression investigation on stability of mxene materials. Computational Materials Science, 196: 10578, 2021. doi: https://doi.org: http://doo.or/1O.... 1016/j.
Nommatsci.2021.11055
Samuel Holt. Zhaozhi Qian, and Mihaela der Schaar. Deep generative symbolic regression. In The Eleventh International Conference on Learning Representations, OpenReview.net, 2023. U

---

Published as a conference paper at ICLR 2025
Alexa-Ale Kamienny, Stéphane d' d'Ascoli, Guillaume Lample, and François Charton. End- to-end symbolic regression with transformers. In Advances in Neural Information Processing Systems, 2022. URL
nips. http:cccppppeeeeeeeeeeeeeeeeeeee
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei. Scaling laws for neural language models. arXiv preprint arXiv:2001.0836/ 2020.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 14/2.1422.. 2014.
William G. 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, 2016. doi:
William G. La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabrício Olivetti de França, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason H. Moore. Contem- porary symbolic regression methods and their relative performance. In Proceed- ings of the Neural Information Processing Systems Track on Datasets and Bench-
marks, 2021. URL otw.::/wwe//aaatttoooooooooooooa www.uugss..eerii...o.. paper fww.s/ pales//pppr/22211
Guillaume Lample and François Charton. Deep learning for symbolic mathematics. arXiv preprint 1911.09912.044111 2019.
Mikel Landajuela, Chak Shing Lee, Jiachen Yang. Ruben Glatt, Claudio P Santi- ago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. A unified framework for deep symbolic regression. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Infor- mation Processing Systems, volume 35, pp. wwtppppppfellle/pppeeeeeeeeeeeee / 2022. URL httpps:/ http:///prreeednnn.......i...... cc
Curran Associates, Inc.,
fil male/dbcabae htte://bbbbefbboocceoossoodeede2 Sat---e88---a-err- Cewerr-Coeeeeeeneeccm pd.f
Juho Lee, Yoonho Lee, Jungtaek Kim, Adam Kosiorek, Seungjin Choi, and Yee Whye Teh. Set transformer: A framework for attention-basse sermutation-invarian neural networks. In Interna- PMLR, 2019.
tional conference on machine learning, pp.
Wenqiang Li, Weijun Li, Lina Yu, Min Wu, Linjun Sun, Jingyi Liu, Yanjie Li, Shu Wei, Yusong Deng, and Meilan Hao. neural-guided dynamic symbolic network for exploring mathematical expressions from data. International Conference Machine Learning, 2024.
Yanjie Li, Weijun Li, Lina Yu, Min Wu, Jinyi Liu, Wenqiang Li, Meilan Hao, Shu Wei, and Yusong Deng. Metasymnet: A dynamic symbolic regression network capable of evolving into arbitrary formulations. arXiv preprint vrXiv:2311.07326 2023.
Z. Li and H. A. Scheraga. Monte carlo-minimization approach to the multiple-minima problem in protein folding. Proceedings of the National Academy of Sciences, 84(19):661 166115 1987. doi: 10. 10.1073/ppas 19.666611
Jiacheng Liu and Siqi Guo. Symbolic regression in financial economics. In The First Tiny Papers Track at ICLR 2023, Tiny Papers @ ICLR 2023, Kigali, Rwanda, May 5, 2023. OpenReview.net. 2023.
Jingyi Liu, Weijun Li, Lina Yu, Min Wu, Linjun Sun, Wengiang Li, and Yanjie Li. Snr: Symbolic bassdr rectifiable learning framework for symbolic regression. Neural Networks, 165:1021-1034, 2023. ISSN 0893-6080. doi: 2023.06.046. URL https: wwwww www.eceeeeeeooooooooooo rect ccmm/sccienne//ttiicl///p
S0893608023003566
Nour Makke, Mohammad Amin Sadeghi, and Sanjay Chawla. Symbolic regression for interpretable scientific discovery. In Big Sie-Data-Antee Analyiins in Astronomy, Science. and Engineering, pp. 26-40. Springer International Publishing, 2022.
12

---

Published as a conference paper at ICLR 2025
Georg Martius and Christoph H. Lampert. Extrapolation and learning equations. In 5th International Conference on Learning Representations, Workshop Track Proceedings. OpenReview.net, 2017.
URL http:///pppnneeeeeeeeeeeeeeeeeee
net.
Yoshitomo Matsubara, Naoya Chiba, Ryo Igarashi, and Yoshitaka Ushiku. Rethinking symbolic regression datasets and benchmarks for scientific discovery. Journal of Data-centric Machine Learning Research, 2024. BRR.httsss//opeeeeiiiet net. wtt.::/eeee.ommom.ooooroooorr/
Trent McConaghy. FFX: Fast, Scalable, Deterministic Symbolic Regression Technology, pp. 235- 260. Springer New York, New York, NY, 2011. doi: 10.1007/978-1-6611--77---
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Ed- ward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):108711092 1953. doi: 10.1063/1 1.006//..6991144 1699114.
T. Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Cláudio P. Santiago, Daniel M. Faissol, and Brenden K. Petersen. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. In Advances in Neural Information Processing Systems, volume 34, pp- 24912-24923, 2021.
URL http:sss//oocedding..eeeeee http::/www....pspppppppeeeeeeeee cc/paper/2021
Jorge
Nocedal and Stephen J. Wright. Springer Series in Operations Research and Springer New York, 2 edition, 2006. doi:
Numerical Financial
Optimization.
Engineering.
Zuzana Oplatkova and Ivan Zelinka. Symbolic regression and evolutionary computation in setting an optimal trajectory for a robot. In 18th International Workshop on Database and Expert Systems Applications (DEXA 2007), pp. (68-172, 2007. doi: 10.1 1,.11099EXX.2000 2.XX....0..5.
Brenden K. Petersen, Mikel Landajuela, T. Nathan Mundhenk, Cláudio Prata Santiago, Sookyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. 9th International Conference on Learning Repre- sentations, ICLR 2021. OpenReview.net, 2021. URL/ttpss:/ wwwp://iee//.r.eie.......... net forum? idd=msQshOBBQG
Markus Quade, Markus Abel, Kamran Shafi, Robert Niven, and Bernd R. Noack. Prediction of dynamical systems by symbolic regression. Phys. Rev. E, 94:012214, 2016. doi: 10.1103/
PhysRevE.94.012214,
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 4439-4447. PMLR, 2018. URL http://
proceedings prw.s v80 hat00/aahoolaa...m
Michael D. Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data.
Science, 324(5923):81 1-885, 2009. doi: 10. 10.1126/science. 165893.
Michael D. Schmidt and Hod Lipson. Age-fitness pareto optimization. In Genetic and Evolutionary Computation Conference, GECCO 2010, Proceedings, pp. 543-544. ACM, 2010. doi: 10.1145/ 1830483.1830584
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas. Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE, 104(1): 148-175, 2015.
Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discovering governing equations via monte carlo tree search. arXiv preprint 2rX.v:2005. 13134, 2022.
Silviu-Marian Udrescu and Max Tegmark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 2666),, 2020. doi: doi: 1000100.1266sciadvvay223.3
Nguyen Quang Uy Nguyen Xuan Hoai, Michael O O'Neill, Robert McKay, and Edgar Galván López. Semantically- crossover in genetic programming: application to real-valued sym- bolic regression. Genet. Program. Evolvable Mach, 12(2)911-19, 2011. doi: 10.1007/ $10710-010-9121-2 RRI.ihtps::///e :ttp::/:://oi..o//1...00.00 1000..1077/000
13

---

Published as a conference paper at ICLR 2025
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nel- son, Eric Jones, Robert Kern, Eric Larson, C Carey, Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, António H. Ribeiro, Fabian Pedregosa, Paul van Mul- bregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-277 2020. doi: 10.1038/st 1992-199-2------ 022---011---2222
David J. Wales and Jonathan P. K. Doye. Global optimization by basin-hopping and the lowest en- ergy structures of -ennard-jones clusters containing up to 110 atoms. Journal of Physical Chem- istry A, 101:5111-511 1997. doi: https://doi.org/ 10. 102111j9770994..
Yiqun Wang, Nicholas Wagner, and James M. Rondinelli. Symbolic regression in materials science. MRS Communications, 930,933-0, 2019. doi: 10. 10.1557/mrc22019.85
Herbert S Wilf. generatingfiunctioolloyy CRC press, 2005.
Graham Woan. The Cambridge handbook of formulas. Cambridge University Press, 2000.
Matthew Wormington, Charles Panaccione, Kevin M Matney, and D Keith Bowen. Characterization of structures from x-ray scattering data using genetic algorithms. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357 (1761):2827-2844)) 1999.
Yang Xiang, DY Sun, W Fan, and XG Gong. Generalized simulated annealing algorithm and its application to the thomson model. Physics Letters A, 233 1997.
14

---

Published as a conference paper at ICLR 2025
A IMPLEMENTATION DETAILS OF PARFAM
In this section, some further implementation details are discussed.
A.
REGULARIZZTION OF THE DENOMINATOR
Since aim for simple function representations, i... for sparse solutions 0 RT. the regular- ization term R()) is of great importance. If we parameterize a rational function Q :: ::R. + in one dimension by
i== a,at
Q(x) (z)
(6)
Ei= b,r*
with a R='+1 and b E IRd++, the following problem occurs: Since for any 7 R \ {0} and (a, b) E RI x Rd it holds that Q(a,)) (z) = Q (2) (z), the parameters canmot be uniquely determined. Although the -on-uniquensss of the solution is not problem in itself, it shows that this parameterization is not the most efficient, and, more importantly, the regularization will be bypassed since 7 be chosen arbitrarily small. We address this issue by normalizing the coefficients of the denominator, i.e., we use b &
Teo rather than b. In other words, instead of defining rational
functions by equation 6 we consider
Σe' 0,æ0
((a,)(x))
= beot
Note that using the 2-norm and not the 1-norm is important since we regularize the coefficients using the 1-norm. To illustrate this, let b
|| 1-6 b.
Case p == When p 1, we have |ر|| for any b E Rdd This demonstrates that i is not regularized anymore and, consequently, also b is not regularized. In essence, this choice of does not promote sparsity in the solution.
ll ف ll1. This expression favors sparse solutions, as encourages the elements of i to close to zero, thus promoting regularization and sparsity in the solution.
Case = 2: In contrast, when p = 2. we have |18|
A.2 MISCELLANBOUS
In general, we look for rational functions Qi whose numerator and denominator polynomials have degree greater than in order model functions like of exp(2x2). However, for some base functions, such as exp, / sin, cos, higher powers introduce redundancy, for instance, exp(r2)² = exp(2r2). To keep the dimension of the parameter space as small as possible without limiting the expressivity of ParFam, we allow the user to specify the highest allowed power of each chosen base function separately. In our experiments, we set it to for all used basis functions: exp, cos and V.
In addition, to ensure that the functions generated during the optimization process are always well- defined and we do not run into an overflow, we employ various strategies:
To ensure that /Q(z)) (z) well-defined, i..., Q(x) 2 0 for x the data set, we use Q(x)] instead.
To avoid the overflow that may be caused by the exponential function, we substitute by the approximation min miinfexp(((r)) exp(10) + [Q(x))), which keeps the interesting regime but does not run into numerical issues for big values of Q(r)). However, adding ((x)) ensures that the gradient still points to a smaller Q(x).
To stabilize the division and avoid the division by 0 completely, we substitute the denomi- nator by 10-3 if its absolute value is smaller than 10-5
Implementing further base functions can be handled similar way as for the square root if they only defined subset of TR. or are prone cause numerical problems.
15

---

Published as a conference paper at ICLR 2025
Symbolic recovery rate
0.8
Accuracy solution rate
0.6
Basin-hopping
0.4
0.6
BFGS
Differential evolution
0.2
0.4
Dual annealing
L-BFGS
10
10
102
Training time s
Training
Symbolic solution rate
(b) Accuracy solution rate
Figure 6: Symbolic solution and accuracy solution rate (percentage of data sets with R > 0.999 for the test set) of ParFam with different optimizers on the subset of the Feynman problems displayed Table
B OPTIMIZER COMPARISON
As discussed in the main paper, ParFam needs to be coupled with a (global) optimizer to approximate the desired function. This section compares different global optimizers, underpinning our decision to use basin-hopping. We tested the following optimizers, covering different global optimizers and local optimizers combined with multi-start:
L-BFGS with multi-start (Nocedal Wright 2006
BFGS with multi-start (Nocedal & Wright 2006
Basin-hopping (Wales & Doye
1997)
Dual annealing (Xiang et al.
1997)
Differential evolution Wormington et al. 1999)
We conducted the experiments on a random subset of 15 Feynman problems, which are listed in Table in Appendix F For each of the 15 problems, we ran ParFam with each optimizer for seven different random seeds and different numbers of iterations. As we solely compare the influence of different optimizers in this experiment, each algorithm. Hence, we are only learning the parameters of one parametric family (So (fo))gR have to since the perfect model parameters result in a parametric family with too many parameters to be optimized in reasonable time and, thus, wasting unreasonable resources. The results are presented Figure 6 These show the superiority of basin-hopping and BFGS with multi-start compared all the other algorithms. While basin-hopping and BFGS with multi-start perform similarly well, it is notable that basin-hopping is less sensitive to the training time (and hence the number of iterations). Therefore, we chose basin-hopping in the main paper, although using BFGS with multi start would have led to similar results.
assume full knowledge of the perfect model parameters
instead of iterating through multiple ones as in the experiments in Section 3 Therefore, omit the problem Feynman-test-
16

---

Published as a conference paper at ICLR 2025
C DATASET CREATION FOR DL-PARFAM
The synthetic data set for the training of DL.PParFam is sampled in the following way. For each dimension 1
I of and each model parameter of ParFam (i.e., maximal degree of the polynomials, base functions, etc.), we sample the same amount of functions from the parametric family fo by sampling 0 IR.. the coefficients however, is not directly sampled from specified distribution, as we have to prevent the sampling of extremely complicated and unrealistic formulas. Therefore, is restricted to sparse, by limiting the number of coefficients allowed be non-zen zero. Specifically, for each polynomial involved, number between and 3 is chosen, which determines the number of coefficients of the denominator and numerator polynomials that are allowed to be non-zero.
Choosing random subsets of the parameters to non-zero, however, might produce function that is covered by more restrictive set of model parameters. For example, if the degree of the numerator of the output rational function Qx+1 is supposed to be 2, but the coefficients of all monomials of degree 2 are with degree for the numerator of Qx+1- This hinders the training of the network since a data set have multiple correct labels.
0, then the function is also covered by the "smaller" set of model parameters,
For this reason, we restrict the random subset of coefficients to include all coefficients that are necessary such that the generated function can be modeled by smaller parametric family. This process ensures that each function has as a target the "smallest" model parameter necessary to describe them and, therefore, for each input the network there is unique target.
Note that, it might be possible that some parts of the function can be simplified due to mathematical cquivalences. However, many of these equivalences depend highly on the basis functions used, so concentrating on these would mean to overfit to specific functions which we aim to avoid.
After choosing the non-zero entries of 0, sample them from (O, 9) and sample I1, TN U(1, 5). The next step is to compute y/s = fo (z;). As it is possible that the sampled function fe is well-defined on the whole domain, some of the values y, might be NaN oc. Furthermore, we want to restrict ourselves to functions with yı E [-1000, 1000) to ensure that the sampled functions are reasonable. Therefore, it necessary to filter these functions afterward.
Our procedure for filtering is follows: For 6 times we resample all x; for which |y:| > 1000 or y/ was NaN. This helps to keep functions that only have a singularity somewhere in the domain but otherwise interesting. If after times, there are still points in the domain with [yx| > 1000 y; NaN, we remove the generated function fe from the data set. In Table we present a random subset of the generated expressions.
17

---

Published as a conference paper at ICLR 2025
Table 2: Example formulas of the synthetically generated expressions for the pre-training of DL- ParFam
Formula
1.0( 111s2 2.612
1.0( 1.333x 2.5110+2.12
-2.665.88
+ 4.032mme
1.00 1.001.80.2gg-00 -0.674)
O.502z 0.887 0.462
3.841
243
0.874
0.4882
4.074g0=1
0.65120=2
2.420
0.340 44,000
3.646.z
4202 00.744g
2.226 cos
0.753rz + 20.3. 0375mg
0.5..1222
0.3922g cos
(1.657a
59 -0..3331 =
cos
0
0.29720 Xn
0.669x
-2 80820 - 2.206
0.5
-0.743
2 5.888 - 22208
B.5
0.506x2 0.819
0.505
0.819 10.108
7885m
0.8 11 o cos
0.176ku
0.575 cos
0.1762g
0.100z +0002240..7000 40.714gg
son 10.18
0.109
40.7
40.7742g
-0.000
3.487
0.666x1222
1.626e
-0.082g
0.3821
0.780zg
0.483
1.0(S.598x2
1.0
4.9652g
1.02
1.762x, cos
33
1.48219446
7.804a 2 1.411egmn 68dngiig -0 178e 2,300 2 817egiy 050 2 68.0
2.
12.12 U.0snag
00.225
0.931mg +0.36627
D IMPLEMENTATION DETAILS OF DL-PARFAM
In this section, we give further details on DL-ParFam.
Embedding The input (æt, U:):= may vary in the sequence length N and the dimension of I; E R.. To deal with a varying sequence length, we apply a Set Transformer (Lee et al. (019)) since the ordering of the input-output pairs is important. To deal with the varying dimension of 24. < m < 9, we embed them in R Since we noticed that linear embedding improves the performance only marginally, we opted for a simple embedding by appending zeros.
Target encoding The targets of the network are the model parameters for ParFam. We focus here on five different ones: The degrees the numerator and the denominator of first layer rationals Q1...., Qk. the degrees of the numerator and the denominator of the second layer rational Qu+1 and the functions 91...., 9k used. We opted for one-hot encoding for all five parameters.
At inference time At inference we apply the pre-trained network predict probabilities for the different model parameters for ParFam. We then extract m sets of model parameters with the high- est probability. for ParFam to try. Furthermore, we add the model parameters for polynomial and rational functions, since these involve fewer parameters and can be checked quickly. In this work, we set m 3 to speed up the computations. Other works (Kamienny al. 2022) propose to sample from the predicted probabilities to increase the diversity, which we did not consider in this work but might be a useful extension in the future.
18

---

Published as a conference paper at ICLR 2025
Table 3: The training and data parameters for the pre-training of DL-ParFan
Maximal Degree First Layer Numerator
2
Maximal Degree First Layer Denominator
2
Maximal Degree Second Layer Numerator
4
Maximal Degree Second Layer Denominator
3
Base functions
v cos, exp
Data parameters
Maximal potence of any variable
3
(i.e., 1 is excluded but 2 is allowed)
Number data pairs (Iz; y1) per function
200
Minimal dimension of s;;
1
Maximal dimension of x;
9
Optimizer
ADAM
Step size
0.0001
Batchsize
1024
Training parameters
Gradient clipping
1
Hidden dimension
256
Number inducing points (Set Transformer)
128
Number heads per multi-head
4
Number layers encoder
8
Number layers classifier
4
In Table we show the generalization gap of the SET Transformer trained for DL-ParFam, which arises by training it synthetic data and then applying it to the Feynman data sets. We evaluated how often the SET Transformer correctly predicts model parameters for synthetic training datasets and the Feynman dataset, considering the top & most likely predictions.
Top
Top 3
Top 5
Top 10
Synthetic
31.4%
50.2%
61 61.8%
71.2%
Feynman
30.4%
38.0%
40.5%
45.6%
Table Percentage of data sets for which the SET Transformer of DL--ParFam predicted the correct model parameters on the synthetic and the Feynman datasets.
The results indicate that while the SET Transformer used for DL-ParFam generalizes well to OOD data (the Feynman data sets) for its top predictions, there is room to optimize the synthetic train- ing data further to improve its generalization. Note that "correct model parameters" refer to those spanning the parametric family with the minimal number of parameters covering the target function. Thus, DL-ParFam can sometimes recover the correct function without using the exact "correct" model parameters.
E EXPRESSIVITY OF PARFAM
Our goal in this section is to compute an estimate for the ratio q//d, To this end, we first consider the number binary trees of complexity I, which denote by by. Note that there exist n binary trees with 0 internal nodes, thus, by n. A binary tree with internal nodes can be created by combining two binary trees, whose internal nodes sum up to 2- with binary operator as the root. This results in the following recurrent formula:
- bi, by -1 for 1. 0. bo n.
(8)
19

---

Published as a conference paper at ICLR 2025
ParFam represent any unary-binary tree where each path from the root leaf has most one unary node. Thus, it holds co = n. A new tree with / internal nodes can be created by either adding unary node as the root of binary tree with / - 1 internal nodes or by adding binary node the root of two trees that can be represented by ParFam with - internal nodes together. Therefore,
a = kb + b b 911111 for / > 0, n.
Lastly, we analyze the number of all unary-binate trees with internal nodes, including those not representable by ParFam. There exist n unary-binary trees with 0 internal nodes, so dg = n. A new unary-binary tree with internal nodes can be created by either adding unary node the root of unary-binary tree with 7-1 internal nodes or by adding a binary node as the root of two unary-binary trees with /- internal nodes together:
di = + bb, di, 11-11-1, for dg = n.
(10)
Given specific values for k. b. and n, for the first few values of I, shown Figure 7 in Appendix E.I To compute the ratio q in general, however,
the formulas for b1, ct, and d, to calculate /
need compute an explicit exact or approximate formula for e and dj. We start by deriving an approximate formula for q. Therefore, we compute the generating function (Wilf, 2005) of bt, following the ideas Lample Charton 2019). For the ease of notation we set di 0 for all / < O. -h-4hmn
Lemma E.1. For the generating function of (bi)eez given by B(2) Eiez B(2)
2bz
The proof can be found in Appendix
E.2
Based on this formula, we can determine the generating function C(x) of [e)) The proof follows similar argumentation as the one for (xx) and can be found in Appendix E.3
Lemma E.2. For the generating function of (ci)ee given by C(z) = E C(s) = B(s)k= 2bx
To derive an approximation of & we compute the singularity with the smallest absolute value of C,
see Theorem 5.3.1 in Wilf 2005). The singularities of are the zeroes of
p(z) (bz(k=B(2)) + n) 2k + 2kz 4bzn 4bnz.
(11)
Computing zero of p equivalent finding a solution of
(1- 2kz - 4bnz)² = 4k2 2 - 4bnz).
(12)
The solutions of this equation can be computed using Cardano's formula. Therefore, we define
= n kn2 10bk2 34.3
36bb-n
(13)
:= -2* -96b>km
2
14066 - -
864b4--y
and
862n² + 13bkn + 16
V 23+0++R R2 +2
(14)
276486245m3
With the quantities
S= v R+ VV and T= VR-VV.
(15)
the roots of p given by
bn k 31.2
31 S
S
bn + k
iv3
42
2 (S T)
(16)
2
3k2
S + T
bn + k:
iv3
T3
(S T).
2
3k2
2
20

---

Published as a conference paper at ICLR 2025
Evaluating p on 21, 22, and Ty shows that 21 and 22 the zeros of p, if choose the branches of the - in S and 7 such that Ve' 6 Choosing different branches results in the same zeros of but a different numbering. In Appendix E.4 we show that |z|| < |122 and, thus, |z]| is the singularity with the smallest absolute value of C. As discussed before, this allows us compute an explicit approximation of ct. Theorem E.3. It holds that
10 (VAT(i 4 +13
3 AV
+ 0 (17)
C1
2bxr 1+1
8/4T(1 1)5
U1
1)
for some constants vo, V1 IR (depending on b, and n).
In Appendix we additionally compute the exact formulas for uo and 01. Moreover, we show the approximations and the true values e in Figure Sa revealing that already the first-order approxima- tion of q is quite close to the exact one.
Since interested in et also need approximation of d. We again start by computing the di
generating function:
Lemma E.4. For the generating function of (di)lez. given by D(2) = Σ, diz',
D(2) 1-ks /k2 2 4bm)a+1
2b
The proof can be found in Appendix E.6 Lample & Charton (2019) also proved this lemma, aiming for an approximation of di afterwards well. However, their calculation included some small typos. For this, we derive a different approximation in Appendix E.7 and consider a higher order approximation, which closer as shown Figure 8b
+ 2bn + 2 B2n2
Theorem E.5. For 11,
holds that
k2
di
12
+
3r
+ +0 (18)
2br
(In(++)3
(
8/
E.1 VISUALIZATION OF THE RATIO c//4
0.8
G/ld,
0.4
0=
7
10
20
25
30
Complexity f
Figure 7: Ratio of and d, for b 4 and k 3 for different values of n and 0..,, 30.
computed using equation and equation TO
E.2 PROOF LEMMA
Proof. We start with multiplying equation g with zt and then summing over all / E Z {0) which yields for the left hand side
Σ b12
biz
bo B(2)
(19)
1200
122
21

---

Published as a conference paper at ICLR 2025
For the right hand side, we get
(I
by, bi -
z = =b
bi, bi -
et = Σ bu-
bis
IGE
===0
(20)
bz Σ b1-12'-1
Σ bizt
bzB(z)²
Therefore,
bzB(2)² - B(z) - = 0,
(21)
which is solved by
B1,2(2) 1 v1 4bbm
(22)
2hs
Since know that B(0) bo and lim=4o B1(z) and lim=-0 B2(z) n, the generating function is given by
B(z) Bz(2) - v1 4ban
(23)
2bz
E.3 PROOF OF LEMMA E.2
Proof. As before, we start with multiplying the recurrence relation for o equation 9 and then sum over all I E Z {0). For the left-hand side this yields again
Σ cas' 190
Σ
East cu C(z) -n.
(24)
For the right-hand side, get
h Σ
b
1
+
Σ
Σ
a 1h---111 d
140
1x0
V1=0
kz
b 0
+ +b Σ
Σ
Cr an-1111) 2
120
11=0
(25)
Σ Ea-id
kzz(z) +b
EacIsti) a
kzz((z) + bz
kzB(2)) + ba((2)2
Together this yields
b=C(z)² C(2) + kz2(22 + =0,
(26)
which can be solved by
C12(22) 1 v1 4bz +n)
(27)
2bz
As before, we know that C(0) CD = n. Therefore, C # C₁ since lim Cz(z) = oo and, thus,
C(=) C2(z) 1 4bz (kz (2)) + n)
(28)
2hz
22

---

Published as a conference paper at ICLR 2025
E.4 PROOF THAT [21] < [
The definition of #1 and 22 is given in equation T6
We want prove that æ1 is the singularity of C
with the smallest absolute value (except for = 0). (We already proved in the main paper that
and z2 are the only relevant singularities.)
Proof. Remember that we chose the branches of the -oots in S and T such that Veto e/ Start with observing that R, V E ce and, therefore, R + vV = xe for some À > 0 and ④ E (T/22, .)). Thus, S = e//3 and é/3 E (Λ/6, π a/3), so Re(S), Im(s) > and /m(s) > ( Next observe that I' is the complex conjugate of S and, therefore, S+T 2Re(5) and S - T= 12Im(5). All together this yields
S
bn + t
iv3
|zz|
(5 -
2
3k2
2
bn + k
iv3
Re(S)
(i2Im(S))
3k2
bm + k 3k2
Re(S)
+ v3lm(S)
> Re(S) + bn + k
+
v (a)(
(29)
3k:2
Re(S) + bn + k
1
+
v3
Re(s)
3k2
v3
2Re(S) + bn + k
3k22
bn + k
12Re(8)
|21|
3k2
E.5 PROOF OF THEOREM 2.1
Proof. Our plan is use Theorem 5.3.1
Wilf 20055 stating that
m
(1-3---1) +0(5
B
Σ Uys
0(s 00-m-1-----))
(30)
where [2']f(2) denotes the 1-th coefficient (corresponding 2²) of the series of powers of
and u=2)= = Ej=a 1 2/s)² is analytic on a disk [z] < |s| + r7 for some 17 > 0 and B & N and
m E N. We define
R(2) = 2k2 + 2k2v1 4bam 4bnz
(31)
such that
C(2) 1- R(=)
(32)
2bz
Hence, the coefficients of the power series of R will allow to compute the coefficients of the power series of C. To compute the power series of R we set
222) (1
1-12
(33)
27
Since R only has the two singularities 2; and x2 and the product of two functions which are analytic point an is itself analytic 20 know that has at most the singularities æn and 22. To
23

---

Published as a conference paper at ICLR 2025
that is analytic around Tı, calculate
(()) =
- 1/2 2kz + 2k V1 4bzn 4bnz
(34)
11
V (1 -2k2 4bnz + 2k 4bzm) (1 2k 4bnz
2k / 4b2n)
(1
2 3-112
I1
V(1 2k 4bnz 22 4bzm)
(35)
2 1/2 V (1 2k 4bmz)² - (2k:2)2(1 4ben)
(36)
I1
2k: 4bnz -k / 4bzn)
The numerator in the last step is the same polynomial arising in equation I and, thus, has the roots Fi, 22, and Tz. Therefore, we can factorize it and continue the computation to
p(l -1- -/)
o(z) (1
_) 2 31
1/2
(37)
(1
2kiz Abmz 2k v1 4b2n)
(u(1 xz
2,)(1 2.)
(38)
(1 2k 4bnz 2k v 4ben)
Since (1 2k 4bnz 2k /1 4bzn) has that v is analytic in and, thus, 22 is the only singularity of v. Furthermore, since we have shown in Appendix [E.4 that [x|| [122 we know that there is some r7 >> 0 such that analytic on some disk |z| < |mn] + Thus, we can use equation t0 get
with μ
21, this shows
[=)) u027 (1-1/--1) + 112 eier"( 3/2 (1-31((1 + O1ET'1-712)
(39)
First we need to determine vo and U1 which compute by developing around æ1;
v(2) + 21) v(x1) 210' 2/21).
(40)
So, we compute
21 )1 22
= v(x1) / (1 2kx - 4bnix - 22 vI 4bx 17)
22
(41)
and
Vj =- =x1)
(42)
2x1
2 //111 22 -11)(1 - 2)[[ 2kz1 4bnr, 2k v 4bzz
D223
222
23
(43)
v 111 21 (1
21 22
+ 01 2(1 2k | 4bnxt
(44)
2k:0 /1
4dx 3/22
4bnkx
(-2k: 4bm 2h: 4bxrn +
(45)
v1
4bx1
Furthermore, we can use the formula (Wilf 2005)
a 1)
1 1-
a(a + 1)
21 +
0
(46)
24

---

Published as a conference paper at ICLR 2025
10f3
G
First order approximation approximation
10
+
dr
First order approximation
www.soo---------- ------1----11111
Second order approxtmation
10a4
10a
Number of trees
10
10
10
10a
1034
1010
104
10
5
15
25
30
10
15 Complexity
25
Complexity /
(a) The values of and its approximations puted using equation 55 (first order approximation) and equation 5] (second order approximation) for 6 = 4, k
(b) The true values of c and its approximations computed using equation 668 (first order approxima- tion) and equation 69 (second order approximation) for b 4. l 3, and n 5.
3. and n 5.
Figure 8: The true values of a and of de and their approximations computed in Theorem 2.1 and 2.2
for b 4, h = 3, and = 5.
to compute
1/2 1)
13//2
1/2- (1/2
3
[1
+ 0(()))
7/2
[(-1/2)
21
1
8
(47)
and
1 (1-3/2-1) = l(-3/2)
- 3/2
1-5/2
[1 +(
+ +0
(48)
TV
where we used that - - and I ((-3/2)
1/ 3
Using the approximation from equation 47 and equation
48
yields
[=1)(=)) = x =
1
3
+
3
0
3
8
mv
(49)
1
= =- 00
+
13
We can now further compute, for >>
0,
q = (=(00(((=)))>
R(=)
1
((()
1
-11
x
1
552))
2b
2
2h
2b AT + 1)3
(50)
[1(2)
(51)
26
1
x1 1-1l00
vol
+
3
1- 1-1
3
0
117/2
25
/ + 1)3
8 8 + 1)3
4 T( 1)3
(52)
E.6 PROOF OF LEMMA E.4
Proof. As for by and & we get for the left-hand side
Σ
d
d,2' do D(2)
(53)
100
le2
25

---

Published as a conference paper at ICLR 2025
and for the right-hand side
k: Σ
d +6
di
Σ
Σ di d o
140
140
kz
di dl-12' + E
Σ
di, d )
at
\1=0
(54)
kz Σ 21-12'-11 +
Σ di-12")
[dde))
1zD(2)) + bz Σ Ediieott)
Σ
Ediz))
kzD(2) + bzD2
Together this yields
bzD(z)² + (kz (z) n
(55)
which is solved by
D1,((2)= -kk - 1) + v(kz -1)2
4bbn
(56)
2bz
As before, using D(0) = n yields
D(2) D2(x) -(ki 1) v k2 (2k: + 4bm)z
(57)
2bz
E.7
PROOF THEOREM 2.2
As for Theorem 2.1 we plan to use Theorem 5.3.1 in Wilf 2005) which yields that
Σ Uys 1 --3 1-1111 )) + 01---11----3) (1
2/5)
(58)
1=0
where [='][(=) denotes the 1-th coefficient (corresponding to 2²) of the series of powers of f 8o and v(2) = Ei= 0y 2/ 2 is analytic a disk [a] |s| + 17 for some 17 > 0 and B & N and
m E N. Similar to the proof of Theorem 2.1 in Appendix E.5 we start by defining
R(=) = /2 (2k + 4bm)z
(59)
such that
D(2) 1- k- R(2)
(60)
2b
To be able to define first have compute the singularities of R with the smallest absolute value. The singularities of R are the zeros of p(z) = 2²22 (2k + 4bn)z + which are
2k 4bn + (2k + 4bn)2 4k2 + 2bn + V(k: 2bn)2 K
112
(61)
21
22
Since h + 2bn > O and (k + 2bm)2 2 2.2 > 0, we know that |ral > |r for all b, k, n > 0. Thus, we define
(()) = (1 - X
(62)
r2
r1
where À = kk//rr which is analytic on disk | < |r2l + n for some v > since Iral > |rel-
Choosing m in equation 58 only need to compute uo and e which can do by developing v around r2: m(z) u(r2) (((r2((= r2) v(r))= ravl (rz)( =/r2)
(63)
26

---

Published as a conference paper at ICLR 2025
which shows that
20 = v(r2) = AV1 r2/r1
(64)
and
Xr2
U1 -r2u((())
(65)
2/1 r1
since
A
(2)
(66)
- r 2
Plugging Up and by and the approximations for the binomial coefficients from equation 47 and equation 48 into equation 58 yields
[=)) = r XV1 r2/ri(- -rr/nili
3
+ -
dr
+
10
V4Tl3
2 2/11 1n Tl
1 +
= = A r2/m1
V4T13
(67)
We can now further compute, for ( >> 1,
di d=)=)(()))
- 1
|59 R(z)
1
|++
R(x)) ~ ~:2""=1-000 a r2/mi
1
(68)
2b
2h
2 / (1 + 1)3
or
di
1
(2(t))
(69)
2h
2h
r2
V4+(++1)
8445(1+1}3
r2
111 Ar
2 1-
12 m
4
-(+11}
(70)
which perfectly approximates dj.
27

---

Published as a conference paper at ICLR 2025
F SRBENCH DATASETS
SRBench La Cava et al. 2021) comprises two types of problem sets: the ground-truth datasets, which include the Feynman datasets (Udrescu & Tegmark 2020) and the Strogatz datasets (La Cava et al. 2016), and the black-box datasets. The ground-truth datasets are synthetically generated and follow known analytical formulas, while the black-box datasets consist of real-world data, where the existence of concise analytic formulas is uncertain, and they are often high-dimensional. In the following, we describe the different data sets in more detail.
FEYNMAN DATA SET
The Feynman data set consists of 119 physical formulas taken from the Feynman lectures and other seminal physics books (Udrescu & Tegmark 2020). Some examples can be found in Appendix FF The formulas depend a maximum of 9 independent variables and composed of the elemen- tary functions +, *, exp, log, sin, cos, tanh, arcsin and arccos. Following La Cava et al. (021), we omit three formulas containing arcsin and arccos and one data set where the ground- truth formula missing. Additionally, since the data sets contain more data points than required for the reconstruction of the equations and this abundance of data slows down the optimizer, we only consider a subset 500, for the experiments without noise, and 1,000, for the experiments with noise, data points of the training data for each problem. We only use the full data sets for EndTTEEn Kamicnny et al. 2022 to perform the bagging described in their paper.
Table 5 shows a random subset of the Feynman data set. The complete Feynman data set can be downloaded from the Penn Machine Learning Benchmarks (MIT license).
28

---

Published as a conference paper at ICLR 2025
Table 5: Random subset of 15 equations of the Feynman problem set (Udrescu & Tegmark 2020
Name
Formula
hw
Feynman-I1l-4-33
2T (exp
hws Tk
2T Ent
Feynman-111-8554
= sin? 2
h
1pyyman--l--1--4
V
Bmom cos (0)
Feynman-1l-24-17
π²
d
2x BJ zgmom
Feynman-11-34-29h
h
Feynman-I-12-5
Efgz
mi71 + m222
Feynman-I-18-4
y
mi + 122
ch?
Feynman-1-38-1
kmqq
Tkba V
Feynman-l-39-22
gmz
Feynman-I-400
exp
Tkb
Feynman-II44 1-43-311
y Tkbmob
Feynman-I-8-14
V (-x1 + x2)² + (-y1 + y2)
Gmn
Feynman-I-99 1-9-18
(-I1 + x2)² + (-81 + V22 + (-21 22)²
az
122
Feynman-testt
2m
3 H3 + 72
Feynman-test-I 18
y
8TG
29

---

Published as a conference paper at ICLR 2025
F.2 STROGATZ DATA SET
The Strogatz data set introduced by La Cava et al. (2016) is the second ground-truth problem set included in SRBench (La Cava 2021 It consists of 14 non-linear differential equations
describing seven chaotic dynamic systems in two dimensions, listed in Appendix F.2 set contains 400 samples.
Each problem
Shows the complete Strogatz data set. It can be downloaded from the Penn Machine Learning Benchmarks (MIT license).
Table 6: The Strogatz ODE problem set (La Cava et al. 2016).
Name
Formula
0.5x* x: + 20
Bacterial Respiration
ý = 0.5.
an + 10
④ sin (x)) + 0.5 sin (Ir - y)
Bar Magnets
sin (y) 0.5 sin (æ: y)
0.05... sin (v)
Glider
= x (y)
④ x²
2xy 3x
Lotka- Voltera interspecics dynamics
xy
y² 2y
x
- x+1
y +
Predator Prey
0.075y
cos
(x)) cot (1)
Shear Flow
(0. sin² (v) + cos2 (v)) sin (x)
=
10a"
10w
+ 10y
van der Pol oscillator
F.3 BLACK-BOL BOX DATA SETS
The black-box data sets comprise 133 complicated, real-world data sets for which the underlying formula there even is one, is unknown. The problems are often high-dimensional. As mentioned in the limitations, ParFam is not yet suited for extremely nigh-dimensional datasets, so we restrict our focus to black-box datasets with dimensionality of 10 or fewer, which leaves us with 77 out of the total 122 datasets.
G BENCHMARKED ALGORITHMS
For experiments the Feynman (Udrescu & Tegmark 2020) and the Strogatz (La Cava ct al. 2016) datasets we decided to besidessssid the algorithms benchmarked La Cava al. 2021 algorithms which showed state-of-the-a art performance in other experiments and pre- training based methods which uploaded their code and weights to have good comparisons for Par- Fam and DL-ParFam.
PySR PySR (Cranmer 2023) uses a multi multi-ppopultion evolutionary algorithm, which consists of unique wwwpssiiippllyy-piiiiiii loop its search algorithm. Furthermore, it involves noise
30

---

Published as a conference paper at ICLR 2025
filter to become more robust and applicable to real data. In the experiments, used 6 populations of size 50 and 500 cycles per iteration. The formulas are restricted to have a maximum depth of 10 and a maximum size of 50. These are the default parameters shown in the tutorial, which also performed the best on a grid hyperparameter search on the first 20 problems of the Feynman dataset. Our experiments use the PySR package version 0.18.1 (Apache license 2.0).
uDSR uDSR (Landajuela et al., 2022) is a combination of the recursive problem simplification introduced in Udrescu & Tegmark (2020) DSR (Petersen al. 2021). genetic programming, pre- training, and linear models. It can, therefore, be seen as the extension of the combination of DSR and genetie programming introduced in Mundhenk et al. (021). The experiments Landajuela et al. (2022) performed in their paper on the Feynman and Strogatz were limited to 24CPUh and 2,000,000 evaluations. For this reason, we reran uDSR our machines with the limits specified by La Cava let al. (202T): 8CPUh and 1,000,000. Because of this, the performance we report is slightly worse. We performed the same hyperparameter search as explained for PySR and picked, in the end, their default hyperparameters (specified in the config regression.gp.jee commit: 2069d4e). Our experi- ments use the official DSO github repository BSD 3-Clause License).
End EndToEEd (Kamienny et al.) 2022) utilizes a pre-trained transformer to predict the symbolic form of the formula directly with estimations of the constants well. To the best of our knowledge, this is the only pre-trained method that is able to handle 9 dimensions and arbitrary constants-which is necessary for the Feynman dataset-for which the model weights available. Instead of the hyperparameters specified in the paper (100 bagging bags, 10 expression trees per bag, and 10 trees that are refined) we use 500 bagging bags, 10 expression trees per bag, and 100 trees that are refined since we experienced a better performance with these. Our experiments use the official github repository (Apache 2.0 license).
NeSymRes NeSymRes (Biggio et al., 2021 works similarly to EndToEnd (Kamienny et al. 2022). however, they only predict the skeleton of the expression, leaving placeholders for all the constants. These are learned afterward using BFGS. Since NeSymRes only works for dimensions to 3, we run it on restricted subset of the Feynman and Strogatz datasets, the results are shown in Appendix N Since we experienced the best performance with them, we performed the experiments with the hyperparameters setting specified in Biggio et al. (202T): beam size of 32 and restarts. Our experiments use the official github repository (MIT license).
H MODEL PARAMETER SEARCH FOR PARFAM
The success of ParFam depends strongly on good choice of the model parameters: The set of base functions 91..... Th and the degrees di and df, i {l...., k 1}. of the numerator and denominator polynomials of Q1, .... Qk+1, respectively. On the one hand, choosing the degrees very small the set of base functions narrow might restrict the expressivity of ParFam too strongly and exclude the target function from its search space. On the other hand, choosing the degrees too high or a very broad set of base functions can yield a search space that is too high-dimensional to be efficiently handled by global optimization method. This might prevent ParFam from identifying even very simple functions.
To balance this tradeoff. we allow ParFam iterate through many different choices for the hy- perparameters describing the model. The user specifies upper bounds on the degrees df and of of the polynomials and the set of base functions 91. 9k- ParFam then automatically traverses through different settings, starting from simple polynomials to rational functions to more complex structures involving the base functions and ascending degrees of the polynomials. The exact pro- cedure is shown in Algorithm [ Note that we refer to the rational functions QI...., Qu. which will be the inputs to the base functions, as the 'input rationals' and, therefore, describe the degrees of their polynomials by DeglnputNummrrtor and DeglnputDenommiatoo.. Similarly, we denote the degrees of the polynomials of the output rational function Qx+1 by 'DegOuputNumerrtto and DDegOutputDenominato
31

---

Published as a conference paper at ICLR 2025
Algorithm 1: Traversal of the model parameters
Input:
Maximal Degree Input Numerator d maxxin"
Maximal Degree Output Numerator d maxxoout
Maximal Degree Input Denominator de maxxinn
Maximal Degree Output Denominator
de max
Maximal number of base functions 6max
Set base functions Gmax = 191.....gk]
Output: List of model parameters L. that define the models ParFam iterate through.
Let C= be an empty list.
// Start with polynomial model:
D ("DeglaputNumerrtoor" O, "ee@uuputNumerator" dmax,out> 'DeglapuiDenominttor: 0, "DegOuuputDennoiiaaoo: "baseFunctions": [}} L.append(Do)
Dp =
// Continue with purely rational models with different degrees:
for daut d2 lmax.ouu do di du do D. "neglnnuuNuumrrtooo::: 0, 'DegOutpuNummerator:: dbut. 'DeglnputDenooinatorr:
for d
0, DeggOupputDenominator:: deut. 'baseFunctions': []] L.))
end
end
// Include different combinations of base functions:
10
for b= 1 to max do
for Le for d'
d2 mux do to d maax, do
11
for di = 0 to de max do
for d = 1 to d max.in do
for B as a list with b elements of Gmax do
/ / Note that base functions can contained in B multiple times.
D = I'DegInputNumerrtoor: d,, 'DegOutpuNNumerator', d
DeglnpuuDenominatoo's di, 'DegOutputDenominator. d2
basstinntooss': B} Caapeed(()) D))
end
19
end
end 20
21
end
22
end
end
24 return C
This strategy is comparable the one proposed by Bartlett et al. (202)), called "Exhaustive Sym- bolic Regression". There, they iterate through a list of parameterized functions and use BFGS to identify the parameters. To create the list of parametrized functions, they construct every possible function using a given set of base operations and a predefined complexity. Notably, this results in more than 100,000 functions to evaluate for one-dimensional data, with the same of base func- tions as we do, but without cos. Our algorithm, however, only needs to search for the parameters of around 500 functions since covers many at the same time by employing the global optimization strategy.
Due this high complexity, Bartlett dimensional problems and, thus, could benchmark their algorithm only on Feynman-1-6-2a (y cxp(8? /2) /2pi). the only one-dimensional problem from the Feynman data set (Udrescu & Tegmark 2020). This example shows the benefit of employing global search in the parameter space:
2023 state that they merely concentrate one-
32

---

Published as a conference paper at ICLR 2025
While ParFam needs five minutes of CPU time to compute the correct function, Bartlett et al. (2023) need 33 hours (150 hours, if the set of possible functions is not pre-generated).
HYPERPARAMETER SETTINGS SRBENCH GROUND-TRUTH PROBLEMS
The hyperparameter settings for the SRBench ground-truth problems are summarized Table 7
Table The model and optimization parameters for the SRBench ground-trui problems for ParFam and DL-ParFam
Maximal Degree First Layer Numerator
2
Maximal Degree First Layer Denominator
2
Maximal Degree Second Layer Numerator
4
Model parameters
Maximal Degree Second Layer Denominator
3
Base functions
V
cOs, exp
Maximal potency of any variable
3
(i.e., 21 is excluded but uï is allowed)
Global optimizer
Basin-hopping
Local optimizer
BFGS
Optimization parameters
Maximal number of iterations global optimizer
10 (1 for DL-ParFam)
Maximal data length
1000
Regularization parameter
0.001
J ADDITIONAL PLOTS FOR THE SRBENCH GROUND-TRUTH RESULTS
Feynman
Strogatz
ParFam DL-Par
X
uDSR
AlFeynman
AFP FE
DSR
gpleam
GP-GOMEA
ITEA
EPLEX
Operon
SBP-GP
BSR
Target Noise
EndToEnd
0.0
FFX
MRGP
0.1
25
75
50
Solution Rate (%6)
Solution Rate (%)
Figure 9: Symbolic solution rate on both SRBench ground-truth data sets separated.
33

---

Published as a conference paper at ICLR 2025
Feynman
Strogatz
ParFam DL-ParFam
PySR
uDSR
AlFeynman
AFP FE
DSR
AFP
gpleam
GP-GOMEA
ITEA
EPLEX
Operon
SBP-GP
BSR
Target Noise
EndTeEnd
0.0
FFX
MRGP
0.1
0.5
1.0
0.0
1.0
Accuracy Solution
Accuracy Solution
Figure 10: Accuracy solution rate (percentage of data sets with R 0.999 for the test set) on the
SRBench ground-truth problems separately.
SRBench ground-truth problems
R? Test
Model Size
Training Time (s)
ParFam
DL-ParFam
x
AlReynman
AFP
DSR
AP
Target Noise
spleam
GP GOMEA
0.01
ITEA
EPLEX
0.1
SBP-GP
BSR
EndToEEd
FFX
MRGP
ex
1-10-1
1-10-
1- 10
1-10-12 1
1-10-1
10:
10a
102
10a
10a
10a
10a
10*
Figure 11: Results on the SRBench ground-truth problems. Points indicate the median test set performance all problems. The R² Test for AlFeynman is missing the plot since SRBench used a higher precision data type, such that AlFeynman achieved a median R² greater than 10 16
K SENSITIVITY ANALYSIS FOR X
In Table K we present the results for ParFam on the ground-truth SRBench data sets for different values of X. Note, that this has been done afterwards as a sensitivity analysis and not to choose 0.001 was based on theoretical considerations and prior observations
X. Our selection of on toy examples while debugging ParFam. Table & shows that ParFam is robust with respect to X. Surprisingly, the complexity of the learned formulas increases for increasing This counterintuitive phenomenon might be due to various reasons. First, the norm does not enforce sparsity but favors it, since is only a proxy for it. So, lower norm does not necessarily imply a lower sparsity. Furthermore, the enumeration through the model parameters breaks the monotonous influence of the regularization. For example, a smaller parametric family might have been the best for lower regularization parameter.
34

---

Published as a conference paper at ICLR 2025
Table Results of ParFam on the ground-truth SRBench data sets for different values of A.
X
Accuracy solution rate
Symbolic solution rate
Complexity
0,0001
94.7%
50%
227
0.001
91.7%
55.6%
131
0.01
94.7%
52.2%
243
L BLACK-BOL -BX DATA SETS: DL-PARFAM
The current implementation of DL-ParFam can only handle independent variables at most. For this reason, benchmark DL-ParFam against the other algorithms on the black-box data sets, which have at most 9 independent variables. The result on the remaining 50 data sets can be seen in Figure 12
R Test
Model Size
Treining Time (s)
*Py5R
*pariam
-Operon
*GI-GOMEA
*D
RidgeRidge
*FEAT
ATTEA
P
Adallcost
RandomForest
MLP
*FFX
LGBM
Linear *MRGP
FBSR
*All eynman
-0.25 0.00 0.25
1.00
101
10a
107
Figure 12: Median R², formula complexity, and training time on the 50 black-box problems from SRBench (La Cava et al. 202I) with at most independent variables. The asterisk indicates that is symbolic regression method.
M COMPARING PARFAM TO EQL ON SRBENCH
As described in the introduction, EQL Martius & Lampert 2017) is the closest method to ParFam, since both make use of non-linear parametric models to translate SR continuous optimization problem. Because of this similarity, we believe that it is important to also show numerical compar- isons between these two. Even though Sahoo et al. (2018) extended EQL include the division operator, also their version of EQL (EQL with division) is not able to express the square root, log- arithm, and exponential, which is why we created reduced version of the ground-truth SRBench, which omits all equations using any of these base functions, In total, this covers 96 formulas. The results of EQL with division on these can be seen in Table
To ensure a fair comparison for EQL, we first tried to run it using the default learning parameters and model parameter search recommended by the authors. However, since EQL. will then quickly use up the computing budget given by SRBench (8 hours of CPU time) we tested EQL for multiple dif- ferent hyperparameters on the first 20 problems from SRBench. We then chose the pest-performing hyperparameters (epoch factor: 1000, penalty every: 50, maximal number of layers: 4, 11: 10) and reran the whole benchmark. This, together with the initial run using the recommended parameters, gives two formulas per equation. In the results shown in Table 9 we chose the formula with the better R² on the validation data set. Note, that we did not make use of the information, that the square root, logarithm, and exponential never parts of the formulas when running ParFam, i.e.,
35

---

Published as a conference paper at ICLR 2025
included these base functions in the dictionary. In experiments used the code from the official EQL github repository (GNU General Public License v3.0).
Table Results of ParFam and EQL with division (Martius & Lampert, 2017 Sahoo et al. 2018 on the 96 SRBench ground-truth equations, which do not include the square root, logarithm, and exponential.
Accuracy solution (R² 0.999)
Symbolic solution
ParFam
93.8%
69.8%
EQL
75%
16.7%
N COMPARING DL-PARFAM TO NESYMRES ON SRBENCH
Since NeSymRes (Biggio et al. 2021) allows most independent variables, we compare it with Da-ParFam a corresponding subset of SRBench. We perform the experiments with the same settings as the experiments shown Section the same hyperparameters for DI-PPaFFam, which reported in Appendix [ The hyperparameters for NeSymRes results are summarized in Table I0] which shows that DL-ParFam outperforms NeSymRes in the symbolic solution rate well as the accuracy solution rate while being more than 20 times faster, even though the model was trained for dimensions to 9.
reported in Appendix G The
O COMPARISON WITH BAYESIAN OPTIMIZATION
The goal DL-ParFam is to simplify the model-parameter search, which usually done by per- forming grid search for ParFam, where we start with testing simple configurations first and slowly increase the complexity of the parametric families. Another standard approach to accelerate the hyper-parameter optimization is Bayesian optimization (Shahriari et al. (015). To evaluate the im- pact by pre-training a SET Transformer first, to guide the model-parameter selection, we compare the performance of DLmPPrFamm with ParFam with Bayesian optimization using Gaussian processes to guide the model-paramete selection and ParFam with structured grid search. The Bayesian op- timization searches through the same model parameters ParFam with grid-search (Algorithm [[]. ie., the values shown in Table 7 Table I presents the results on the ground-truth problems of SRBench. While Bayesian hyperparameter optimization manages to speed up the training as well, DL-ParFam outperforms it with respect to symbolic solution rate and training time.
P NGUYEN BENCHMARK
To compare ParFam with SPL (Sun et al. 2022) and NGGP (Mundhenk et al. 2021). which are the
current state-of-the-art some SR benchmarks (like Nguyen (Uy al. 2001 but results of them on SRBench were reported, we evaluate ParFam on Nguyen. Interestingly, we observed that the original domain on which the data was sampled is big enough specify the functions, as
Table 10: Results on the subset of the SRBench ground-trutt problems containing only expression with most 3 variables. Following SRBench terminology, training time refers to the time each algorithm requires to compute a result for specific problem, which corresponds to inference time for pre-trained methods.
Symbolic solution
Accuracy solution (R² > 0.999)
Training time (in s)
DL-ParFam 58.1%
87.1%
53
NeSymRes
48.7%
59.0%
1389
36

---

Published as a conference paper at ICLR 2025
Symbolic solution rate
Accuracy solution rate
Training time
Bayesian (max. 50 calls)
34.9%
85.3%
7678s
Bayesian (max. 500 calls)
38.0%
89.1 %
10937s
DL-ParFam
45.9%
83.5%
234s
Grid search
55.6%
93.2%
12860s
Table 11: Comparison of symbolic solution rate, accuracy solution rate, and training time on the ground-truth problems of SRBench for ParFam with model-parameter selection guided by Bayesian optimization, a SET Transformer, and grid search.
ParFam was able to find simple and near indistinguishable approximations to the data that are not the target formula. For example, it found 0.569.. 0.742 sinn1.241.r 2.059) 1.655 instead of (2² cos(x) 1, since both are almost identical on the domain [-1, 1]. this reason, we extended the data domain for some of the problems. The results for the Nguyen data set can be in Table [2 We used the hyperparameters shown in Table 13 Following (unn et al. (2022). from whom take the results of the competitors,
sin and exp the standard basis functions for ParFam and add / and log for the problems 8, 11, and 8 Note that the formula Nguyen-11 can not be expressed by ParFam and hence the symbolic solution rate is
37

---

Published as a conference paper at ICLR 2025
Table 12: Results on the Nguyen benchmarks. The results for ParFam are averaged over 6 inde- pendent runs. The results from SPL (Sun et al., 2022), NGGP (Mundhenk et al., 2021), and GP (a genetic programming-bann -baeed SR algorithm) are taken from Sun et al. (2022)
Benchmark
Expression
ParFam
SPL
NGGP
GP
N
3.
x +
100%
100%
100%
99%
Nguyen-2
x²
+ x² + 2 +
100%
100%
100%
90%
Nguyen-3
w + x² + x² x2 +
100%
100%
100%
34%
Nguyen-4
x + x + x + x + x.² + x
100%
99%
100%
54%
Nguyen-5
sin (x²2 cos (x) -
83%
95%
80%
12%
6
sin (I) sin (x² + x)
83%
100%
100%
11%
Nguyen-7
log (x + 1) + log (x² + 1)
100%
100%
100%
17%
Nguyen-8
v
100%
100%
100%
100%
Nguyen-9
sin (mo) + sin (z?)
100%
100%
100%
76%
Nguyen-10
sin (mo) cos (x1)
100%
100%
100%
86%
Nguyen-11
x
0%
100%
100%
13%
Nguyen-12
af - z + 0.5z - 21
100%
28%
4%
0%
Nguyen-1
3.39x + 2.12x² + 1.78x
100%
100%
100%
0%
Nguyen-2
0.48a + 3.39. + 2.12. + 1.78
100%
94%
100%
0%
Nguyen-5
sin (x2) cos (x) 0.75
83%
95%
98%
1%
Nguyen-8
V1.23.
100%
100%
100%
56%
Nguyen-9
sin (1.52)) + sin (0.5.²)
100%
96%
90%
0%
Average
91.2%
94.5%
92.4%
38.2%
Table 13: The model and optimization parameters for the Nguyen benchmark.
Maximal Degree Input Numerator
2
Maximal Degree Input Denominator
0
Maximal Degree Output Numerator
6
Model parameters
Maximal Degree Input Denominator
0
Base functions
cos, exp (v. log)
Maximal potence of any variable
6
Global optimizer
Basin-hopping
Local optimizer
BFGS
Optimization parameters
Maximal number of iterations global optimizer
30
Regularization parameter À
0.1
38

---

Published as a conference paper at ICLR 2025
Table 14: Results on SRSD-Feynman easy
Symbolic solution
Accuracy solution (R² > 0.999)
Training time (in
ParFam
76.7%
93.3%
7620
DL-ParFam
60%
66.7%
72
PySR
90%
93.3%
5796
uDSR
66.7%
86.7%
7470
Table 15: Results SRSD-Feynman medium medium
Symbolic solution
Accuracy solution (R²² > 0.999)
Training time (in s)
ParFam
47.5%
80%
8478
DL-ParFam
40%
65%
138
PySR
62.5%
85%
8022
uDSR
40%
82.5%
7740
Q SRSD FEYNMAN
Matsubara et al. 2024) introduced the SRSD-Feynman benchmarks (Creative Commons Attribution 4.0 International) building the Feynman datasets (Udrescu & 'Tgmmark] (2200) to substitute the artificial ranges and coefficients imposed by the original datasets with the physical ones, resulting in the same formulas with coefficients and variables ranging from 10 10-30 to 10 The best-performing algorithms in their benchmark are PySR (Cranmer 2023) and uDSR (Landajuela et al. 2022), which why we ParFam and DL-ParFam against those two.
We ran the experiments for the easy and medium difficulty for ParFam and DL-ParFam and the two best-performing algorithms from Matsubara et al. (2024) since these seem to the strongest the field currently, also supported by our experiments on SRBench. We use the same settings as for the Feynman data sets (SCPUh and 1, 1,000,000 function evaluations) and the same hyperparameters for all algorithms as shown in Appendix and G
The results are shown in Table 14 and 15 Note that our results differ from those reported in Mat- subara et al. (2024) since we normalize the data for each algorithm. The results show that ParFam performs slightly worse than PySR but better than uDSR. Furthermore, DL-ParFam is only slightly worse than its competitors, while being up to 50 to 100 times faster. It is particularly interesting to see that Da-PPaFamm has reasonable performance even on data sampled from a very different domain than the one it trained on, probably due to the normalization of the input data during its training.
39


Reviews:
[{"review_id": "IY9V4DxyvQ", "reviewer": "Reviewer_JBKY", "summary": "The authors propose a SR method that fixes a rational function structure and optimizes it using MC-based optimization methods. They also present a version of this technique that uses pre-trained transformers as a starting point. Results are presented on many SR benchmark problems and include an analysis of expressivity.", "questions": "> Typical formulas from the Feynman dataset Udrescu & Tegmark (2020), for instance, have a complexity of around 10.\n\nWhat does this mean?\n\nTable 1: It wasn't clear to me why parfarm's ability to cover possible expressions would increase as $n$ increased. Is it just because $n$ artificially inflates the proportion of copies of identical function forms with swapped leaves in $d$? Or something else?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and to my knowledge a fairly novel approach to SR that is well-contextualized. The authors have made very, very extensive experimental comparisons, although many of these are sent to the appendix. The authors also do a good job of providing a slightly deeper theoretical justification for their work than in typical SR papers which I enjoyed reading. I think the paper makes a good contribution to the field in showing another avenue for SR that revisits established global optimizers with a unique functional structure that can still find ground truth solutions with good fidelity.", "weaknesses": "The main weakness of Feynman and Strogatz dataset comparisons is that many of the equation forms are very simple (especially feynman). So, one way to do well on them is to restrict the complexity of models during optimization (or invent a method that happens to search over simple forms). That's why it's important to also consider benchmarks on real world data. The authors are aware of this, but focus the main body of their text on these synthetic datasets and send real-world results to the appendix. I would be in favor of the real-world dataset comparisons from SRBench being more of a main result, especially if they are given an extra page in the revisions. \n\nI would have also liked to see a stronger connection made between the expressive restrictions of their methods (e.g., inability to model deep unary functions) and the distribution of those types of expressions contained in those benchmarks. \n\n- i would have liked to see more motivation for the chosen representation of equations as rational functions earlier on. Why are shallow rational functions a good _hypothetical_ choice for representing any possible expression tree, before introducing your approach to measuring expressivity. \n\n- typically theorems are self contained. So in Theorem 2.1 it would be more clear to define $c$, $x$, $l$, etc. rather than having to find them throughout the text. \n\n- I found S 2.2 on the expressivity of ParFarm interesting, but I had trouble extracting an intuition for how expressivity scales from the description and from Table 1. It looks like $x$ is being reused for two different variables, and it also appears that the authors are only considering expressivity for $\\ell = 1$ i.e. equations with a single function. I would have liked to know how expressivity scales as functions grow in complexity ($l$) which seems equally or more important than scaling by the number of leaves ($n$) and unary functions ($k$).\n\n- DL-ParFam seems to be extremely sensitive to small amounts of noise in the SRBench ground truth problems. I would have liked to see this mentioned in the results/discussion.\n\n- misc\n  - awkward phrasing: \"to avoid the regularization to be circumvented\"", "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": "2maujaS2Xs", "reviewer": "Reviewer_ZXN3", "summary": "In this paper, we propose a symbolic regression algorithm ParFam, which treats symbolic regression problem as a global optimization problem. Compared with the traditional methods that treat symbolic regression problem as a combinatorial optimization problem, ParFAM improves the search efficiency and has the potential to solve high-dimensional problems", "questions": "##### Questions\n\n1. It is not clear to me how the training data for the DL-ParFam is collected. Please describe it in more detail in the article and provide more details.\n2. To reorganize what is the main innovation of this paper, I don't think using the global optimization problem as a knot symbolic regression problem is an innovation of this paper.\n3. The article says that its ability to deal with high-dimensional symbolic regression is stronger than the existing algorithms, so does the article test and compare the ability of each algorithm to deal with high-dimensional symbolic regression problems?\n4. Why does DL-ParFam not compare inference time with pre-trained symbolic regression methods represented by **Neural Symbolic Regression that Scales** and **End-to-end Symbolic Regression with Transformers**?", "limitations": "", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "This paper proposes DL-ParFam, which uses a pre-trained model to predict ParFam's parameters, effectively improving its training efficiency. The authors claim that the training speed is 100 times faster. But the thinking of DL - ParFam and this article is a bit like (https://doi.org/10.48550/arXiv.2309.13705),", "weaknesses": "##### Weaknesses\n\n1. Mentioned in the article with global optimization to solve the problem of symbolic regression is more efficient, but this is not the first to use global optimization algorithm to solve the problem of symbolic regression algorithm, for example, EQL (https://doi.org/10.1109/TNNLS.2020.3017010), MetaSymNet (https://doi.org/10.48550/arXiv.2311.07326). And this paper does not compare with these two algorithms. Please analyze the advantages of ParFam over the above two methods.\n\n2. At the beginning of the **introduction**, the paper mentioned that the simplicity of the expression is very important, but the paper did not evaluate the expression complexity of the algorithm.\n3. Why did the author only add the noise level of 0.01 in the anti-noise experiment? However, in many other symbolic regression algorithms, the anti-noise experiment is more adequate. I think the noise level should be increased to the order of 0.1 to better test the anti-noise ability of ParFam.", "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": "JcbWQJrxnI", "reviewer": "Reviewer_3uKW", "summary": "This paper presents ParFam, a novel approach to symbolic regression that reformulates the discrete optimization problem into a continuous one using parametric families of functions. The DL-ParFam extension introduces a neural-guided approach that incorporates a pre-trained Set Transformer to accelerate the optimization process. The method is evaluated on standard symbolic regression benchmarks and shows competitive performance.", "questions": "- For DL-ParFam, how robust is the pre-trained model to out-of-distribution data? Are there certain types of functions where it consistently fails to provide useful guidance?\n- Could the authors elaborate on how ParFam's performance scales with the number of variables?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The authors provide a thorough theoretical analysis of ParFam's expressivity, demonstrating that it can still represent a substantial proportion of functions of a given complexity. \n- Results on standard benchmark datasets demonstrate ParFam's strong performance in terms of symbolic recovery rate and accuracy.", "weaknesses": "- The pre-training phase for DL-ParFam is computationally expensive, though it only needs to be done once.\n- In the figure 4, the label \"Symbolic Solution Rate\" may be incorrect and should be \"Symbolic Recovery Rate\".\n- As shown in Figure 4, symbolic solution rate drops significantly when noise is introduced. This suggests that the methods are sensitive to even small levels of noise, which could limit their robustness in real-world scenarios where data is often noisy.\n- The method still faces challenges when dealing with high-dimensional problems (more than 10 variables). As the number of parameters increases exponentially with the number of variables, optimization becomes more computationally intensive.\n- The performance of the method is highly dependent on model parameter choices, such as the degree of polynomials and basis functions, requiring some prior knowledge or extensive experimentation to determine the optimal parameters.", "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": "0UXHXpzrSD", "reviewer": "Reviewer_tj4N", "summary": "In this paper, the authors propose a novel architecture, ParFam, to solve symbolic regression through continuous optimization. The authors also introduce the use of a neural network to predict optimal hyperparameters for ParFam.", "questions": "Here are some questions that need to be addressed:\n1. The idea of using a division operator in EQL was proposed in ICML 2018 [1]. Please provide a comparison between ParFam and EQL with division operators to ensure a fair comparison.\n2. Recent works have incorporated neural architecture search with EQL [2] [3]. These works appear to be a superset of the proposed method. The architecture proposed here is manually designed within the neural architecture search space. It is unclear what advantage ParFam has over these variants of EQL.\n3. The idea of predicting hyperparameters using a transformer is interesting. In ParFam, the default method for hyperparameter optimization is grid search. However, Bayesian optimization is a more common approach. Please provide a comparison of the speedup achieved by hyperparameter optimization using a pre-trained transformer versus Bayesian optimization.\n4. The training time of ParFam on the black-box datasets from SRBench is not shown. Please provide this information.\n5. For the limitation related to high-dimensional data, it is claimed that \"the number of parameters grows exponentially with the number of variables.\" However, from Figure 1, it is unclear why the number of parameters would grow exponentially with the number of variables. Please clarify.\n\nReferences:\n\n[1]. Sahoo, Subham, Christoph Lampert, and Georg Martius. \"Learning equations for extrapolation and control.\" International Conference on Machine Learning. PMLR, 2018.\n\n[2]. Dong, Junlan, et al. \"Evolving Equation Learner for Symbolic Regression.\" IEEE Transactions on Evolutionary Computation (2024).\n\n[3]. Li, Wenqiang, et al. \"A Neural-Guided Dynamic Symbolic Network for Exploring Mathematical Expressions from Data.\" Forty-first International Conference on Machine Learning.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The strength of this paper is that the proposed ParFam shows high accuracy on the Feynman and Strogatz datasets from SRBench. The idea of using a neural network to predict optimal hyperparameters is interesting.", "weaknesses": "The weakness of ParFam is that the advantage of ParFam over EQL and some recent EQL variants is not clearly demonstrated. The authors claim that EQL uses linear layers, while ParFam uses rational layers. However, EQL is not limited to using unary functions in basis functions, and if a division function were included in EQL, it would closely resemble ParFam.", "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 (Poster)