# ParFam – Symbolic Regression Based on
Continuous Global Optimization

Philipp Scholl
Ludwig-Maximilians-Universität München
Munich
Germany
<EMAIL_ADDRESS>
&Katharina Bieker
Ludwig-Maximilians-Universität München
Munich
Germany
<EMAIL_ADDRESS>
&Hillary Hauger
Ludwig-Maximilians-Universität München
Munich
Germany
<EMAIL_ADDRESS>
&Gitta Kutyniok
Ludwig-Maximilians-Universität München
Munich Center for Machine Learning (MCML)
Munich
Germany
<EMAIL_ADDRESS>

###### 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 quite complicated and require a lot of
hyperparameter tuning and computational resources. In this paper, we present
our new method _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 powerful global optimizer,
this approach results in an effective method to tackle the problem of SR.
Furthermore, it can be easily extended to more advanced algorithms, e.g., by
adding a deep neural network to find good-fitting parametric families. We
prove the performance of ParFam with extensive numerical experiments based on
the common SR benchmark suit SRBench, showing that we achieve state-of-the-art
results. Our code and results can be found at
https://github.com/Philipp238/parfam.

_K_ eywords Symbolic Regression $\cdot$ Machine Learning $\cdot$ Global
Optimization.

##  1 Introduction

_Symbolic regression_ (SR) describes the task of finding a symbolic function
that accurately represents the connection between given input and output data.
At the same time, the function should be as simple as possible to ensure
robustness against noise and interpretability. This is of particular interest
for applications where the aim is to (mathematically) analyze the resulting
function afterward or get further insights into the process to ensure
trustworthiness, for instance, in physical or chemical sciences (Quade et al.,
2016; Angelis et al., 2023; Wang et al., 2019). The range of possible
applications of SR is therefore vast, from predicting the dynamics of
ecosystems (Chen et al., 2019), forecasting the solar power for energy
production (Quade et al., 2016), estimating the development of financial
markets (Liu and Guo, 2023), analyzing the stability of certain materials (He
and Zhang, 2021) to planning optimal trajectories for robots (Oplatkova and
Zelinka, 2007), to name but a few. Moreover, as Angelis et al. (2023) points
out, the number of papers on SR has increased significantly in recent years,
highlighting the relevance and research interest in this area.

SR is a specific regression task in machine learning that aims to find an
accurate model without any assumption by the user related to the specific data
set. Formally, a symbolic function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$
that accurately fits a given data set
$(x_{i},y_{i})_{i=1,\dots,N}\subseteq\mathbb{R}^{n}\times\mathbb{R}$ is
sought, i.e., it should satisfy $y_{i}=f(x_{i})$ for all data points, or in
the case of noise $y_{i}\approx f(x_{i})$ for all $i\in\\{1,\dots,N\\}$.
Since, in this general setting, there are no assumptions on the structure of
possible models, the search space is infinite-dimensional. In practice,
however, it is necessary to specify the model space in some sense, and all
methods rely in one way or another on certain implicit assumptions in the
modeling process. For example, _genetic programming_ (GP) methods, one of the
most common classes of solution algorithms, require the choice of base
functions that can be combined to build the model (Schmidt and Lipson, 2009;
2010). Nevertheless, unlike other regression tasks, SR aims at 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 and traditional
regression tasks typically assume a strong structure of the unknown function
like linearity or polynomial behavior.

To tackle SR problems, the most established methods are based on _genetic
programming_(Augusto and Barbosa, 2000; Schmidt and Lipson, 2009; 2010;
Cranmer, 2023), but nowadays there also exist many solution algorithms that
make use of other machine learning methods, in particular neural networks
(Udrescu and Tegmark, 2020; Martius and Lampert, 2017; Desai and Strachan,
2021; Makke et al., 2022). However, even though there have been many attempts
with complicated procedures to search through the infinite-dimensional space
of functions, many of them show unsatisfactory results when evaluated on
complex benchmarks: La Cava et al. (2021) evaluate 13 state-of-the-art SR
algorithms on the _SRBench_ ground-truth problems: the Feynman (Udrescu and
Tegmark, 2020) 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 differential
equations. Out of the 13 algorithms evaluated by La Cava et al. (2021), all
algorithms find at most $30\%$ of the formulas of each problem set in the
given time of 8 CPU hours, except for AI Feynman (Udrescu and Tegmark, 2020).
AI Feynman, which is based on recursive function simplification inspired by
the structure of the Feynman data set, is able to recover more than 50% of the
Feynman equations but fails for more than 70% for the Strogatz problems. The
rates are even worse for data sets incorporating noise (La Cava et al., 2021;
Cranmer, 2023). In addition to AI Feynman, we are only aware of one other
algorithm, proposed by Holt et al. (2023) after the benchmark by La Cava et
al. (2021), which has demonstrated superior performance on the SRBench ground-
truth data sets while following the SRBench time limit.

In this paper, we introduce the novel algorithm _ParFam_ that addresses SR by
leveraging the inherent structure of physical formulas and, thereby,
translating the discrete optimization problem into a continuous one. This
grants users precise control over the search space and facilitates the
incorporation of gradient-based optimization techniques. More precisely, we
apply _basin-hopping_ , which combines a global random search with a local
search algorithm (Wales and Doye, 1997). Originally, this algorithm was
designed to solve molecular problems and, thus, is suitable for very high-
dimensional landscapes. The details of ParFam are introduced in Section 2.1.
Notably, despite its straightforward nature, ParFam achieves state-of-the-art
performance on the Feynman and Strogatz data set as demonstrated in Section
3.1. Moreover, this structure enables the simple application of pre-trained
NNs to reduce the dimensionality of the search space. This concept is
exemplified by our prototype, _DL-ParFam_ , introduced in Section 2.2. Through
experimentation on a synthetic data set, we demonstrate that DL-ParFam
significantly surpasses ParFam, cf. Section 3.2.

##### Our Contributions

Our key contributions are as follows:

  1. 1.

We introduce ParFam, a new method for SR leveraging the inherent structure of
physical formulas and, thereby, translating the discrete optimization problem
into a continuous one. This results in the following advantages: (1) Enabling
gradient-based optimization techniques; (2) Efficient but simple and user-
friendly setup; (3) State-of-the-art performance on the ground-truth problems
of La Cava et al. (2021), the Feynman and Strogatz data sets.

  2. 2.

Furthermore, we introduce the extension DL-ParFam, which shows how the
structure of ParFam allows for using a pre-trained NN, overcoming the
limitations of previous approaches.

##### Related work

Traditionally, genetic programming have been used for SR to heuristically
search the space of equations given some base functions and operations
(Augusto and Barbosa, 2000; Schmidt and Lipson, 2009; 2010; Cranmer, 2023).
However, due to the accomplishments of neural networks across diverse domains,
numerous researchers aimed to leverage their capabilities within the realm of
SR. Udrescu and Tegmark (2020), for instance, have employed an auxiliary NN to
evaluate data characteristics. In a similar vein, Martius and Lampert (2017),
Sahoo et al. (2018), and Desai and Strachan (2021) used compact NN
architectures with physically meaningful activation functions such as $\sin$
and $\cos$, enabling stochastic gradient descent for the search of symbolic
functions.

The approach by Petersen et al. (2021), on the contrary, relies on
_reinforcement learning_ (RL) to explore the function space, where a policy,
modeled by a recurrent neural network, generates candidate solutions. Mundhenk
et al. (2021) 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 to individual RL and GP
approaches.

Given the simplicity of sampling functions and evaluating them, several
endeavors have emerged to train NNs using synthetic data to predict underlying
functions. Initial simpler approaches were limited by the variability of the
given data set (Biggio et al., 2020) or of the functions (Li et al., 2022).
However, these limitations can be effectively circumvented by more advanced
approaches using the transformer architecture (Biggio et al., 2021; Kamienny
et al., 2022; Holt et al., 2023).

Apart from the algorithms evaluated by La Cava et al. (2021), Deep Generative
Symbolic Regression (DGSR) by Holt et al. (2023) and unified Deep Symbolic
Regression (uDSR) by Landajuela et al. (2022) are the only algorithms—to the
best of our knowledge—which have been evaluated on the whole Feynman data set
and outperformed the state-of-the-art AI Feynman in the symbolic recovery
rate. Notably, DGSR’s success has only been possible by the computationally
expensive pre-training of an encoder-decoder architecture using RL instead of
gradient-based methods to learn invariances of the functions and an additional
finetuning step during inference by performing neural guided priority queue
training (NGPQT) as introduced by Mundhenk et al. (2021). uDSR builds upon the
already well-performing AI Feynman and adds pre-training, genetic programming,
reinforcement learning, and linear regression to it. Unlike for the SRBench
benchmark (La Cava et al., 2021), Landajuela et al. (2022) evaluate their
method with a time limit of 24 instead of 8 hours, which is why we omit uDSR
from our comparisons in Section 3.

ParFam shares conceptual proximity with EQL (Martius and Lampert, 2017; Sahoo
et al., 2018), as both methods assume a structure of general formulas,
effectively translating SR into a continuous optimization problem. However,
while ParFam aims to guarantee a unique parameterization for each function,
EQL exhibits many redundancies that inflate the parameter space. EQL relies on
the local minimizer ADAM (Kingma and 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 straightforward inclusion of operations like
logarithms, roots, and division within unary operators—in contrast to EQL as
priorly noted by Petersen et al. (2021). Similar to DL-Parfam, Liu et al.
(2023) enhanced EQL with a pre-training step. However, this approach still
suffers from the listed structural limitations of EQL.

##  2 Methods

In the following section, we first introduce our new method ParFam, that
exploits a well-suited representation of possible symbolic functions to which
an efficient global optimizer can be applied. Afterward, we discuss the
extension DL-ParFam, which aims to enhance ParFam by utilizing deep learning
to obtain better function representations.

###  2.1 ParFam

The aim of SR is to find a simple and thus interpretable function that
describes the mapping underlying the data $(x_{i},y_{i})_{i=1,\dots,N}$
without many additional assumptions. Typically, a set of base functions, such
as $\\{+,-,^{-1},\exp,\sin,\sqrt{}\\}$, 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. Since most algorithms make no other assumptions on the
function they are looking for, this approach results in a search space that
grows exponentially in the number of base functions, dimensions of $x$, and
depth of the expression trees.

To reduce the complexity of the search space on the one hand and to obtain
more meaningful results on the other hand, some methods apply filters to
prevent the output of unwanted or “unnatural” functions. For instance,
Petersen et al. (2021) prevent their algorithm from creating compositions of
trigonometric functions as $\sin\circ\cos$ since these are rarely encountered
in any scientific domain. Given that the main idea of SR is to gain knowledge
of scientific processes, such structural assumptions appear to be reasonable.
This is also the motivation for restricting the search space in our approach.
Furthermore, we choose the function space such that it can be represented by a
parametric family, and the proper expression can be found by applying a
continuous global optimizer.

####  2.1.1 The Structure of the Parametric Family

The main motivation for ParFam is that most functions appearing in real-world
applications can be represented by functions of certain parametric families.
More precisely, we assume that they can be written in the form

$\begin{gathered}f_{\theta}(x)\coloneqq
Q_{k+1}(x,g_{1}(Q_{1}(x)),g_{2}(Q_{2}(x)),\dots,g_{k}(Q_{k}(x))),\end{gathered}$
(1)

where $Q_{1},...,Q_{k+1}$ are rational functions, $g_{1},...,g_{k}$ are the
unary base functions, which cannot be expressed as rational functions, like
$\sin$, $\exp$, $\sqrt{}$ etc., and $x\in\mathbb{R}^{n}$ is the input vector.
Moreover, $\theta\in\mathbb{R}^{m}$ denotes the coefficients of the individual
polynomials, i.e., of the numerators and denominators of $Q_{1},...,Q_{k+1}$,
which are the learnable parameters of this family of functions. The degrees
$d^{1}_{i}$ and $d^{2}_{i}$, $i\in\\{1,\dots,k+1\\}$, of the numerator and
denominator polynomials of $Q_{1},...,Q_{k+1}$, respectively, and the base
functions $g_{1},...,g_{k}$ are chosen by the user. Depending on the
application, even specialized custom functions can be added to the set of base
functions. This versatility and its simplicity make ParFam a 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.

The parametric family we consider excludes composite functions such as
$\sin\circ\cos$ similarly to Petersen et al. (2021). This is rooted in the
structure of physical formulas we observe, as can be seen, e.g., in the set of
ground-truth problems from SRBench (La Cava et al., 2021), which consists of
129 physically meaningful formulas and only includes one function that does
not follow equation 1:

$\sqrt{(x_{1}^{2}+x_{2}^{2}-2x_{1}x_{2}\cos(\theta_{1}-\theta_{2}))}~{}\text{\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{udrescu2020ai}{\@@citephrase{,
}}{}, I.29.16)}}.$

Furthermore, the “Cambridge Handbook of Physics Formulas" (Woan, 2000)
contains more than 2,000 equations from the major physics topics, among which
only a handful do not follow the structure of equation 1.

####  2.1.2 Optimization

Restricting the search space to functions of the parametric family given by
equation 1 yields the advantage that we can translate the discrete SR problem
into a continuous one, as now the task is to find the parameters of the
rational functions $Q_{1},...,Q_{k+1}$ such that $f_{\theta}$ approximates the
given data $(x_{i},y_{i})_{i=1,\dots,N}$, i.e., we aim to minimize the average
$l_{2}$-distance between $y_{i}$ and $f_{\theta}(x_{i})$. As we aim for
preferably simple functions to derive interpretable and easy-to-analyze
results, a regularization term $R(\theta)$ is added to encourage sparse
parameters. In total, we aim to minimize the loss function

$L(\theta)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-f_{\theta}(x_{i})\right)^{2}+\lambda
R(\theta),$ (2)

where $\lambda>0$ is a hyperparameter to control the weight of the
regularization. Here, we choose $R(\theta)=\|\theta\|_{1}$ as a surrogate for
the number of non-zero parameters, which is known to enforce sparsity in other
areas, e.g., NN training (Bishop, 2006; Goodfellow et al., 2016). In Appendix
A, we discuss how to deal with the regularization of the coefficients of
rational functions in detail.

Although the SR problem is now transformed into a 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. This
is also discussed by Nocedal and Wright (2006) and 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 and Doye
(1997), which combines a local minimizer, e.g., BFGS (Nocedal and Wright,
2006), with a global search technique inspired by Monte-Carlo minimization as
proposed by Li and Scheraga (1987) to cover a larger part of the parameter
space. More precisely, we use the implementation provided by the SciPy library
(Virtanen et al., 2020). Each iteration consists of three steps:

  1. 1.

Random perturbation of the parameters.

  2. 2.

Local minimization, e.g., with the BFGS method.

  3. 3.

Acceptance test based on the function value of the local optimum.

The basic idea of the algorithm is to divide the complex landscape of the loss
function into multiple areas, leading to different optima. These are 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 a 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, a finetuning routine is
initiated. In this process, coefficients that fall below a certain threshold
are set to 0, and the remaining coefficients are optimized using the L-BFGS
method, starting from the previously found parameters. The threshold is
gradually increased from $10^{-5}$ to $10^{-2}$ to encourage further sparsity
in the discovered solutions. This step has been found to be crucial in
enhancing the parameters initially found by basin-hopping.

###  2.2 DL-ParFam

As discussed in the related work section, there have been multiple attempts in
recent years to leverage pre-training for SR, as synthetic data can be easily
generated. Even though modern approaches are able to handle flexible data sets
in high dimensions (Biggio et al., 2021; Kamienny et al., 2022), they fail to
incorporate invariances in the function space, e.g., $x+y$ and $y+x$ are seen
as different functions, as pointed out by Holt et al. (2023). Holt et al.
resolve this by evaluating the generated function to compute the loss and
update the network using RL. This effectively solves the invariance problem,
as can be seen by their state-of-the-art symbolic recovery rate on the Feynman
data set. 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 to be
optimized using suitable algorithms like Policy gradient methods.

Here, we propose our approach DL-ParFam, which aims to combine the best of
both worlds. The idea is to use an NN that, given a data set
$(x_{i},y_{i}=f(x_{i}))_{i=1,\dots,N}$, predicts a sparsity pattern on the
coefficients $\theta$ of the parametric family in equation 1. This sparsity
pattern or mask specifies which parameters should be variable and learned in
the ParFam algorithm and which can be ignored and set to a fixed value of $0$.
The idea of DL-ParFam is visualized in Figure 1. This approach yields
significant improvements compared to ParFam and other pre-training based SR
methods:

  * •

Compared to ParFam: DL-ParFam strongly reduces the dimensions of the
optimization problem considered in ParFam, effectively reducing the
computation time and success rate for any global optimizer.

  * •

Compared to other pre-training based SR methods: DL-ParFam predicts the
structure of the function, which can be directly compared with the ground-
truth and, thereby, avoids the evaluation on the data grid in every training
step and yields an end-to-end differentiable pipeline. In addition, DL-ParFam
adeptly handles function invariances, as we guarantee that each set of
parameters uniquely defines a function via the structure of ParFam.

Figure 1: Schematic illustration of the DL-ParFam method.

The primary purpose of this subsection is to demonstrate the potential of
utilizing ParFam as a means of structuring scientific formulas beyond the
direct optimization presented so far. Our intent is not to present an
implementation of DL-ParFam in this paper that can rival existing deep
learning-based methods on complex benchmarks like the Feynman data set. This
decision is driven by the myriad challenges inherent in benchmarking, such as
high dimensionality, diverse data point distributions, varying numbers of data
points and dimensions, and a plethora of base functions. While these
challenges can be addressed using deep learning techniques, as demonstrated by
Biggio et al. (2021), Kamienny et al. (2022), and Holt et al. (2023), they
require specialized architectures as transformers. For now, we opt for a
straightforward implementation of DL-ParFam to demonstrate its effectiveness
on synthetic data.

The vanilla implementation, which we consider here, uses a simple fully-
connected feedforward neural network
$NN:\mathbb{R}^{N}\rightarrow\mathbb{R}^{m}$ which takes as input the data
$(y_{i})_{i=1,...,N}$ and outputs a mask $c\in[0,1]^{m}$, where $c_{i}$
represents the likelihood that $\theta_{i}$ is needed to represent the sought
symbolic function, i.e., $c_{i}\approx 0$ indicates that the parameter
$\theta_{i}$ is supposed to be 0 and $c_{i}\approx 1$ indicates
$\theta_{i}\neq 0$. To reduce the dimensionality of the NN, we only take the
output data $y$ of the functions as input to the NN. Thus, we implicitly
assume that the input data is sampled on the same grid
$(x_{i})_{i=1,\dots,N}\subseteq\mathbb{R}^{n}$ for all data sets. To ensure
this, we train the NN on synthetically generated data
$(y^{j}=(y_{i}^{j})_{i=1,...,N},c^{j})_{j=1,...,K}$, where $c^{j}\in[0,1]^{m}$
is some mask and $y_{i}^{j}=f_{\theta^{j}}(x_{i})$ is the output of the
corresponding function evaluated on the fixed grid point
$x_{i}\in\mathbb{R}^{n}$. As a loss function with respect to the parameters
$w$ of the NN, we define

$L_{\mathsf{NN}}(w)\coloneqq\sum_{j=1}^{K}\operatorname{BCE}(NN(y^{j}),c^{j}),$
(3)

where $\operatorname{BCE}:[0,1]^{m}\times[0,1]^{m}\rightarrow\mathbb{R}$
denotes the binary-entropy loss, i.e.,

$\begin{gathered}\operatorname{BCE}(\bar{c},c)\coloneqq\frac{1}{K}\sum_{l=1}^{m}-\left(c_{l}\log(\bar{c}_{l})+(1-c_{l})\log(1-\bar{c}_{l})\right).\\\
\end{gathered}$ (4)

Another important difference from previous approaches, not outlined before, is
that DL-ParFam combines the experience gained through pre-training with the
power of the method ParFam, which is highly competitive on its own. In Section
3.2, we show the value of this combination.

##  3 Benchmark

In Section 3.1, we evaluate ParFam on the Feynman (Udrescu and Tegmark, 2020)
and Strogatz (La Cava et al., 2016) data sets and report its performance in
terms of the symbolic recovery rate, the coefficient of determination $R^{2}$,
and the complexity of the derived formula showing that our simple setup
significantly outperforms most existing methods and reaches state-of-the-art
in SR. In Section 3.2, we study DL-ParFam, revealing the vast potential of
adding pre-training to ParFam.

###  3.1 ParFam

After the introduction of the SR benchmark (SRBench) by La Cava et al. (2021),
several researchers, including Mundhenk et al. (2021), Holt et al. (2023),
Kamienny et al. (2022), and Biggio et al. (2021), have reported their findings
using the SRBench’s ground-truth data sets. These data sets are the Feynman
and the Strogatz data set.

##### Feynman data set

The Feynman data set consists of 119 physical formulas taken from the Feynman
lectures and other seminal physics books (Udrescu and Tegmark, 2020). Some
examples can be found in Appendix C. The formulas depend on a maximum of 9
independent variables and are composed of the elementary functions
$+,-,*,/,\sqrt{,}\exp,\log,\sin,\cos,\tanh,\arcsin$ and $\arccos$. Following
La Cava et al. (2021), we omit three formulas containing $\arcsin$ and
$\arccos$ and one data set where the ground-truth formula is missing.
Additionally, since the data sets contain more data points than required for
ParFam and this abundance of data slows down the optimizer, we only consider a
subset of 500, for the experiments without noise, and 1,000, for the
experiments with noise, data points of the training data for each problem.

##### 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 et al., 2021). It
consists of 14 non-linear differential equations describing seven chaotic
dynamic systems in two dimensions, listed in Appendix D. Each data set
contains 400 samples.

##### Metrics

To ensure comparability with the results evaluated on SRBench, we use the same
evaluation metrics as La Cava et al. (2021). First, we report the symbolic
recovery rate, which is the percentage of equations ParFam recovered. Second,
we consider the coefficient of determination

$R^{2}\coloneqq
1-\frac{\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{N}(y_{i}-\bar{y})^{2}},$
(5)

where $\hat{y}_{i}=f_{\theta}(x_{i})$ represents the model’s prediction and
$\bar{y}$ the mean of the output data $y$. The closer $R^{2}$ is to $1$, the
better the model describes the variation in the data. It is a widely used
measure for goodness-of-fit since it is independent of the scale and variation
of the data. Lastly, we report the complexity of our formula based on the
number of mathematical operations following the definition in SRBench. The
original data sets do not include any noise. However, similar to La Cava et
al. (2021), we additionally perform experiments with noise by adding
$\epsilon_{i}\sim
N\left(0,\sigma^{2}\frac{1}{N}\sum_{i=1}^{N}y_{i}^{2}\right)$ to the targets
$y_{i}$, where $\sigma$ denotes the noise level.

##### Hyperparameters

The hyperparameters of ParFam can be divided into two subsets. The first
subset defines the parametric family $(f_{\theta})_{\theta\in\mathbb{R}^{m}}$,
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 a wide range of potential functions. In
this context, we opt for $\sin$, $\exp$, and $\sqrt{}$ as our base functions.
For the “input rational functions” $Q_{1},\ldots,Q_{k}$, we set the degrees of
the numerator and denominator polynomials to 2. For $Q_{k+1}$, we set the
degree of the numerator polynomial to 4 and the denominator polynomial to 3.
This choice results in a parametric family with hundreds of parameters, making
it challenging for global optimization. To address this issue, we iterate
through various smaller parametric families, each contained in this larger
family, see Appendix E for details. The second set of hyperparameters defines
the optimization scheme. Here, we set the regularization parameter to
$\lambda=0.001$, the number of iterations for basin-hopping to $10$, and the
maximal number of BFGS steps for the local search to $100$ times the dimension
of the problem. Our choice of parameters is summarized in Table 4 in Appendix
F.

##### Results

Following La Cava et al. (2021), we allow a maximal training time of $8$ CPU
hours and a maximal number of function evaluations of $1,000,000$. In Figure
2(a), we present the symbolic recovery rate on both data sets together.
ParFam, AI Feynman, and DGSR exhibit exceptional performance, outperforming
all other competitors by a substantial margin (over $25\%$). It is important
to note that AI Feynman performs particularly well on the Feynman data sets
but fails on the Strogatz data set, as shown in Appendix G, indicating that
the algorithm is tailored to the Feynman data set. Since DGSR was not tested
on noisy data and AI Feynman is strongly influenced by noise, ParFam
outperforms all competitors at a noise level of $\sigma=0.01$. Furthermore,
Figure 2(b) shows the accuracy solution, which is the percentage of formulas
for which $R^{2}>0.999$ holds on the test sets. Here, ParFam outperforms all
competitors with and without noise. It is important to note that Holt et al.
(2023) reported a similar metric but with $R^{2}>0.99$ instead of
$R^{2}>0.999$. However, DGSR achieved a value of $90.95\%$ for this less
strict metric, compared to ParFam’s $97.67\%$. Figure 3 reveals that ParFam
achieves a mean $R^{2}$ significantly better than its competitors, albeit with
slightly more complex formulas. Note that the mean, rather than the median, is
shown in this figure since both ParFam and AI Feynman solve over 50% of the
formulas without noise, causing their median performance to simply reflect
this high success rate. The corresponding plot showing the median and
additional results can be found in Appendix G.

(a) Symbolic recovery rate

(b) Accuracy solution

Figure 2: Symbolic recovery and accuracy solution rate (percentage of data
sets with $R^{2}\\!>\\!0.999$ for the test set) on the SRBench ground-truth
problems (Feynman and Strogatz data sets). Figure 3: Results on the SRBench
ground-truth problems. Points indicate the mean over all problems.

###  3.2 DL-ParFam

In this subsection, we aim to demonstrate the potential of DL-ParFam by
conducting experiments on synthetically generated data sets.

##### Synthetic data sets

To evaluate DL-ParFam, we employ synthetic data sets due to its current
prototype status, which limits its ability to handle complex data sets like
the Feynman set. To generate these, we fix the grid $x_{i}=-10+0.1i$,
$i\in\\{1,...,N=200\\}$, and choose one set of model hyperparameters to define
a specific parametric family $(f_{\theta})_{\theta\in\mathbb{R}^{m}}$. Here,
we choose the base functions $\sin$ and $\sqrt{}$, set the degree of all
numerator polynomials to $2$ and of all denominator polynomials to $0$. Then,
we sample $(\theta^{j})_{j\in\\{1,\dots,K\\}}\subset\mathbb{R}^{m}$ following
the scheme described in Appendix H. For each $\theta^{j}$,
$j\in\\{1,\dots,K\\}$, and each $x_{i}$ we evaluate
$f_{\theta^{j}}(x_{i})=y_{i}^{j}$ to obtain $K$ data sets
$((x_{i})_{i=1,...,N},(y_{i}^{j})_{i=1,...,N},\theta^{j})_{j=1,...,K}$. For
our numerical tests, we create two different data sets: The first includes
2,000,000 equations for training the neural network of DL-ParFam, 10,000 for
its validation, and another 10,000 for its testing. The second set consists of
$100$ formulas to compare DL-ParFam with ParFam. Our hyperparameter choices
are summarized in Table 6 in Appendix I.

##### Pre-training of DL-ParFam

We construct the neural network of DL-ParFam as a feedforward NN with one
hidden layer containing 200 neurons. It is trained as described in Section 2.2
using the Adam optimizer with a learning rate of 0.0001 and 20,000 epochs with
500 batches, which takes less than 4h on a TITANRTX GPU. In addition, to
predict a mask $c\in\\{0,1\\}^{m}$ as described in Section 2.2 we set
$c_{k}=0$ if $NN((y_{i})_{i=1,...,N})<0.2$ and $c_{k}=1$ otherwise.

##### Metrics

To evaluate the performance of the NN, we report two metrics: the covering
score and the average successful cover size. The covering score describes the
percentage of formulas for which the mask includes the non-zero parameters,
i.e., $\theta_{k}^{j}\neq 0$ implies $c_{k}^{j}=1$. The average successful
cover size is the mean over the means of $c^{j}$ across all formulas for which
the NN succeeded at identifying the non-zero parameters. Ideally, this value
should be as small as possible, indicating that the mask size is minimized
while still effectively capturing the non-zero parameters. To understand the
influence of the NN in DL-ParFam, we evaluate DL-ParFam and ParFam against
each other on the second synthetic data set and report the symbolic recovery
rate. Here, we assume for both methods that the perfect choice of model
parameters is known, i.e., the same that were used to create the data sets.
This assumption allows us to assess the relative performance of DL-ParFam and
ParFam rather than evaluating their general performance.

##### Results

The NN reaches a covering score of $91.32\%$ on the test data set with an
average successful cover size of $26.62\%$. This indicates that the pre-
training helps to reduce the number of parameters by almost 3/4 in $91.32\%$
of the formulas. The relative performance of ParFam and DL-ParFam is shown in
Table 1, which reveals that DL-ParFam solves consistently $16$-$32\%$ more
equations than ParFam while requiring only approximately half the time.

Table 1: Relative performance and runtime of ParFam and DL-ParFam on the synthetic data set | ParFam | DL-ParFam
---|---|---
# Iterations | Symbolic recovery | Training time | Symbolic recovery | Training time
50 | 32% | 13s | 48% | 7s
100 | 30% | 24s | 46% | 12s
500 | 43% | 116s | 66% | 55s
1000 | 37% | 221s | 69% | 107s

##  4 Discussion and Conclusion

This work introduces ParFam along with its potential extension, DL-ParFam.
Despite its inherent simplicity, ParFam demonstrates remarkable performance,
as shown in Section 3.1. Furthermore, its adaptable structure makes it highly
versatile for specific application scenarios. While DL-ParFam currently only
exists in a prototype form, it already shows the feasibility and potential of
integrating pre-training—a crucial direction in SR as pointed out by Kamienny
et al. (2022); Biggio et al. (2021); Holt et al. (2023)—into the ParFam
framework.

##### Limitations

While the structure of the parametric family of ParFam is undoubtedly its
greatest asset in tackling SR, it can also be considered its most significant
constraint, given that it imposes a tighter constraint on the function space
compared to other methods. However, Figure 2 illustrates, on the one hand,
that several algorithms theoretically capable of identifying specific formulas
do not always achieve this in practice. On the other hand, it demonstrates
that even if ParFam restricts the function space too much, it still manages to
find a formula that approximates the original function with remarkable
accuracy. Another limitation of ParFam is that optimizing high-dimensional
problems ($>$$10$ independent variables) is computationally expensive, given
the exponential growth in the number of parameters with respect to the number
of variables.

##### Future Work

Subsequent efforts will concentrate on advancing both ParFam and DL-ParFam.
With ParFam, several avenues remain unexplored, encompassing diverse forms of
regularization, alternative parametrizations, and the potential incorporation
of custom-tailored optimization techniques. Nonetheless, our primary focus
will be on DL-ParFam, driven by its promising potential, as evidenced by our
experiments. Numerous design choices await exploration, including data
sampling strategies, choice of loss function, architecture selection, and
more. Existing research in this direction will undoubtedly serve as valuable
guidance (Kamienny et al., 2022; Biggio et al., 2021; Holt et al., 2023). We
anticipate that these advancements will facilitate the deployment of even more
expansive parametric families, thereby mitigating the limitations outlined
earlier.

#### Reproducibility statements

In our repository https://github.com/Philipp238/parfam, we include the results
of our experiments and the code and instructions necessary to use our
algorithms and reproduce all experiments shown in Section 3. Furthermore, we
report all settings used in the experiments in the Appendices F and I.

## Acknowledgments

This publication was supported by LMUexcellent, 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 as well as by the
Hightech Agenda Bavaria.

Furthermore, G. Kutyniok was supported in part by the DAAD programme Konrad
Zuse Schools of Excellence in Artificial Intelligence, sponsored by the
Federal Ministry of Education and Research. G. Kutyniok also acknowledges
support from the Munich Center for Machine Learning (MCML) as well as the
German Research Foundation under Grants DFG-SPP-2298, KU 1446/31-1 and KU
1446/32-1 and under Grant DFG-SFB/TR 109 and Project C09.

## References

  * Angelis et al. [2023] D. Angelis, F. Sofos, and T. E. Karakasidis.  Artificial intelligence in physical sciences: Symbolic regression trends and perspectives.  _Archives of Computational Methods in Engineering_ , 30(6):3845–3865, 2023.  doi: 10.1007/s11831-023-09922-z.
  * Augusto and Barbosa [2000] D. A. Augusto and H. J. Barbosa.  Symbolic regression via genetic programming.  In _Proceedings. Vol.1. Sixth Brazilian Symposium on Neural Networks_ , pages 173–178, 2000.  doi: 10.1109/SBRN.2000.889734.
  * Bartlett et al. [2023] D. J. Bartlett, H. Desmond, and P. G. Ferreira.  Exhaustive symbolic regression.  _IEEE Transactions on Evolutionary Computation_ , pages 1–1, 2023.  doi: 10.1109/TEVC.2023.3280250.
  * Biggio et al. [2020] L. Biggio, T. Bendinelli, A. Lucchi, and G. Parascandolo.  A seq2seq approach to symbolic regression.  In _Learning Meets Combinatorial Algorithms at NeurIPS2020_ , 2020.
  * Biggio et al. [2021] L. Biggio, T. Bendinelli, A. Neitz, A. Lucchi, and G. Parascandolo.  Neural symbolic regression that scales.  In M. Meila and T. Zhang, editors, _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_ , pages 936–945. PMLR, 2021.  URL http://proceedings.mlr.press/v139/biggio21a.html.
  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Information Science and Statistics. Springer New York, 2006.
  * Chen et al. [2019] Y. Chen, M. T. Angulo, and Y.-Y. Liu.  Revealing complex ecological dynamics via symbolic regression.  _BioEssays_ , 41(12):1900069, 2019.  doi: 10.1002/bies.201900069.
  * Cranmer [2023] M. Cranmer.  Interpretable machine learning for science with pysr and symbolicregression. jl.  _arXiv preprint arXiv:2305.01582_ , 2023.
  * Desai and Strachan [2021] S. Desai and A. Strachan.  Parsimonious neural networks learn interpretable physical laws.  _Scientific reports_ , 11(1):12761, 2021.
  * Goodfellow et al. [2016] I. J. Goodfellow, Y. Bengio, and A. C. Courville.  _Deep Learning_.  Adaptive computation and machine learning. MIT Press, 2016.  URL http://www.deeplearningbook.org.
  * He and Zhang [2021] M. He and L. Zhang.  Machine learning and symbolic regression investigation on stability of mxene materials.  _Computational Materials Science_ , 196:110578, 2021.  doi: https://doi.org/10.1016/j.commatsci.2021.110578.
  * Holt et al. [2023] S. Holt, Z. Qian, and M. van der Schaar.  Deep generative symbolic regression.  In _The Eleventh International Conference on Learning Representations_. OpenReview.net, 2023.  URL https://openreview.net/pdf?id=o7koEEMA1bR.
  * Kamienny et al. [2022] P. Kamienny, S. d’Ascoli, G. Lample, and F. Charton.  End-to-end symbolic regression with transformers.  In _Advances in Neural Information Processing Systems_ , 2022.  URL http://papers.nips.cc/paper_files/paper/2022.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * La Cava et al. [2016] W. G. La Cava, K. Danai, and L. Spector.  Inference of compact nonlinear dynamic models by epigenetic local search.  _Engineering Applications of Artificial Intelligence_ , 55:292–306, 2016.  doi: 10.1016/j.engappai.2016.07.004.
  * La Cava et al. [2021] W. G. La Cava, P. Orzechowski, B. Burlacu, F. O. de França, M. Virgolin, Y. Jin, M. Kommenda, and J. H. Moore.  Contemporary symbolic regression methods and their relative performance.  In _Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks_ , 2021.  URL https://datasets-benchmarks-proceedings.neurips.cc/paper_files/paper/2021.
  * Landajuela et al. [2022] M. Landajuela, C. S. Lee, J. Yang, R. Glatt, C. P. Santiago, I. Aravena, T. Mundhenk, G. Mulcahy, and B. K. Petersen.  A unified framework for deep symbolic regression.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 33985–33998. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/dbca58f35bddc6e4003b2dd80e42f838-Paper-Conference.pdf.
  * Li et al. [2022] J. Li, Y. Yuan, and H.-B. Shen.  Symbolic expression transformer: A computer vision approach for symbolic regression.  _arXiv preprint arXiv:2205.11798_ , 2022.
  * Li and Scheraga [1987] 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):6611–6615, 1987.  doi: 10.1073/pnas.84.19.6611.
  * Liu and Guo [2023] J. Liu and S. 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.
  * Liu et al. [2023] J. Liu, W. Li, L. Yu, M. Wu, L. Sun, W. Li, and Y. Li.  SNR: Symbolic network-based rectifiable learning framework for symbolic regression.  _Neural Networks_ , 165:1021–1034, 2023.  doi: 10.1016/j.neunet.2023.06.046.  URL https://doi.org/10.1016/j.neunet.2023.06.046.
  * Makke et al. [2022] N. Makke, M. A. Sadeghi, and S. Chawla.  Symbolic regression for interpretable scientific discovery.  In _Big-Data-Analytics in Astronomy, Science, and Engineering_ , pages 26–40. Springer International Publishing, 2022.
  * Martius and Lampert [2017] G. Martius and C. H. Lampert.  Extrapolation and learning equations.  In _5th International Conference on Learning Representations, Workshop Track Proceedings_. OpenReview.net, 2017.  URL https://openreview.net/forum?id=BkgRp0FYe.
  * Metropolis et al. [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.  doi: 10.1063/1.1699114.
  * Mundhenk et al. [2021] T. N. Mundhenk, M. Landajuela, R. Glatt, C. P. Santiago, D. M. Faissol, and B. K. Petersen.  Symbolic regression via deep reinforcement learning enhanced genetic programming seeding.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 24912–24923, 2021.  URL https://proceedings.neurips.cc/paper/2021.
  * Nocedal and Wright [2006] J. Nocedal and S. J. Wright.  _Numerical Optimization_.  Springer Series in Operations Research and Financial Engineering. Springer New York, 2 edition, 2006.  doi: 10.1007/978-0-387-40065-5.
  * Oplatkova and Zelinka [2007] Z. Oplatkova and I. 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)_ , pages 168–172, 2007.  doi: 10.1109/DEXA.2007.58.
  * Petersen et al. [2021] B. K. Petersen, M. Landajuela, T. N. Mundhenk, C. P. Santiago, S. Kim, and J. T. Kim.  Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.  In _9th International Conference on Learning Representations, ICLR 2021_. OpenReview.net, 2021.  URL https://openreview.net/forum?id=m5Qsh0kBQG.
  * Quade et al. [2016] M. Quade, M. Abel, K. Shafi, R. K. Niven, and B. R. Noack.  Prediction of dynamical systems by symbolic regression.  _Phys. Rev. E_ , 94:012214, 2016.  doi: 10.1103/PhysRevE.94.012214.
  * Sahoo et al. [2018] S. Sahoo, C. Lampert, and G. 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_ , pages 4439–4447. PMLR, 2018.  URL http://proceedings.mlr.press/v80/sahoo18a.html.
  * Schmidt and Lipson [2009] M. D. Schmidt and H. Lipson.  Distilling free-form natural laws from experimental data.  _Science_ , 324(5923):81–85, 2009.  doi: 10.1126/science.1165893.
  * Schmidt and Lipson [2010] M. D. Schmidt and H. Lipson.  Age-fitness pareto optimization.  In _Genetic and Evolutionary Computation Conference, GECCO 2010, Proceedings_ , pages 543–544. ACM, 2010.  doi: 10.1145/1830483.1830584.
  * Udrescu and Tegmark [2020] S.-M. Udrescu and M. Tegmark.  AI Feynman: A physics-inspired method for symbolic regression.  _Science Advances_ , 6(16):eaay2631, 2020.  doi: doi:10.1126/sciadv.aay2631.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye.  Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms.  _Journal of Physical Chemistry A_ , 101:5111–5116, 1997.  doi: https://doi.org/10.1021/jp970984n.
  * Wang et al. [2019] Y. Wang, N. Wagner, and J. M. Rondinelli.  Symbolic regression in materials science.  _MRS Communications_ , 9(3):793–805, 2019.  doi: 10.1557/mrc.2019.85.
  * Woan [2000] G. Woan.  _The Cambridge handbook of physics formulas_.  Cambridge University Press, 2000.
  * Wormington et al. [1999] M. Wormington, C. Panaccione, K. M. Matney, and D. K. 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–2848, 1999.
  * Xiang et al. [1997] Y. Xiang, D. Sun, W. Fan, and X. Gong.  Generalized simulated annealing algorithm and its application to the thomson model.  _Physics Letters A_ , 233(3):216–220, 1997.

##  Appendix A Implementation details

In this section, some further implementation details are discussed.

###  A.1 Regularization of the denominator

Since we aim for simple function representations, i.e., for sparse solutions
$\theta~{}\in~{}\mathbb{R}^{m}$, the regularization term $R(\theta)$ is of
great importance. If we parameterize a rational function
$Q:\mathbb{R}\rightarrow\mathbb{R}$ in one dimension by

$Q(x)=Q_{(a,b)}(x)=\frac{\sum_{i=0}^{d^{1}}a_{i}x^{i}}{\sum_{i=0}^{d^{2}}b_{i}x^{i}}$
(6)

with $a\in\mathbb{R}^{d^{1}+1}$ and $b\in\mathbb{R}^{d^{2}+1}$, the following
problem occurs: Since for any $\gamma\in\mathbb{R}\setminus\\{0\\}$ and
$(a,b)\in\mathbb{R}^{d^{1}+1}\times\mathbb{R}^{d^{2}+1}$ it holds that
$Q_{(a,b)}(x)=Q_{(\gamma a,\gamma b)}(x)$, the parameters cannot be uniquely
determined. Although the non-uniqueness of the solution is not a problem in
itself, it shows that this parameterization is not the most efficient, and,
more importantly, the regularization will be bypassed since $\gamma$ can be
chosen arbitrarily small. We address this issue by normalizing the
coefficients of the denominator, i.e., we use $\tilde{b}=\frac{b}{||b||_{2}}$
rather than $b$. In other words, instead of defining rational functions by
equation 6, we consider

$Q_{(a,b)}(x)=\frac{\sum_{i=0}^{d^{1}}a_{i}x^{i}}{\tfrac{1}{||b||_{2}}\sum_{i=0}^{d^{2}}b_{i}x^{i}}.$
(7)

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
$\tilde{b}=\frac{1}{||b||_{p}}b$.

Case $p=1$: When $p=1$, we have $||\tilde{b}||_{1}=1$ for any
$b\in\mathbb{R}^{d^{2}}$. This demonstrates that $\tilde{b}$ is not
regularized anymore and, consequently, also $b$ is not regularized. In
essence, this choice of $p$ does not promote sparsity in the solution.

Case $p=2$: In contrast, when $p=2$, we have
$||\tilde{b}||_{1}=||\frac{b}{||b||_{2}}||_{1}$. This expression favors sparse
solutions, as it encourages the elements of $\tilde{b}$ to be close to zero,
thus promoting regularization and sparsity in the solution.

###  A.2 Miscellaneous

In general, we look for rational functions $Q_{i}$ whose numerator and
denominator polynomials have a degree greater than $1$ in order to model
functions like $x_{1}^{2}\exp(2x_{2})$. However, for some base functions, such
as $\exp,\sqrt{},\sin,\cos$, higher powers introduce redundancy, for instance,
$\exp(x_{2})^{2}=\exp(2x_{2})$. 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 1 for all used basis functions:
$\exp,\cos$ and $\sqrt{}$.

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 $\sqrt{Q(x)}$ is well-defined, i.e., $Q(x)\geq 0$ for all $x$
in the data set, we instead use $\sqrt{|Q(x)|}$.

  * •

To avoid the overflow that may be caused by the exponential function, we
substitute it by the approximation $\min\\{\exp(Q(x)),\exp(10)+|Q(x)|\\}$,
which keeps the interesting regime but does not run into numerical issues for
big values of $Q(x)$. However, adding $|Q(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 denominator by $10^{-5}$ if its absolute value is smaller than
$10^{-5}$.

Implementing further base functions can be handled in a similar way as for the
square root if they are only defined on a subset of $\mathbb{R}$ or are prone
to cause numerical problems.

##  Appendix B Optimizer comparison

As discussed in the main paper, ParFam needs to be coupled with a powerful
(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 and Wright, 2006]

  * •

BFGS with multi-start [Nocedal and Wright, 2006]

  * •

Basin-hopping [Wales and 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 2 in Appendix C. 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, we assume full knowledge of the perfect model
parameters for each algorithm. Hence, we are only learning the parameters
$\theta$ of one parametric family $(f_{\theta})_{\theta\in\mathbb{R}^{m}}$
instead of iterating through multiple ones as in the experiments in Section
3.1. Therefore, we have to omit the problem Feynman-test-17 since the perfect
model parameters result in a parametric family with too many parameters to be
optimized in a reasonable time and, thus, wasting unreasonable resources. The
results are presented in Figure 4. These show the superiority of basin-hopping
and BFGS with multi-start compared to 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.

(a) Symbolic recovery rate

(b) Accuracy solution rate

Figure 4: Symbolic recovery and accuracy solution rate (percentage of data
sets with $R^{2}\\!>\\!0.999$ for the test set) of ParFam with different
optimizers on the subset of the Feynman problems displayed in Table 2.

##  Appendix C Example Feynman problems

Table 2 shows a random subset of the Feynman data set. The complete Feynman
data set can be downloaded from
https://space.mit.edu/home/tegmark/aifeynman.html.

Table 2: Random subset of 15 equations of the Feynman problem set [Udrescu and Tegmark, 2020]. Name | Formula
---|---
Feynman-III-4-33 | $y=\dfrac{h\omega}{2\pi\left(\operatorname{exp}\left({\tfrac{h\omega}{2\pi Tkb}}\right)-1\right)}$
Feynman-III-8-54 | $y=\sin^{2}{\left(\dfrac{2\pi E_{n}t}{h}\right)}$
Feynman-II-15-4 | $y=-Bmom\cos{\left(\theta\right)}$
Feynman-II-24-17 | $y=\sqrt{-\dfrac{\pi^{2}}{d^{2}}+\dfrac{\omega^{2}}{c^{2}}}$
Feynman-II-34-29b | $y=\dfrac{2\pi BJzgmom}{h}$
Feynman-I-12-5 | $y=Efq_{2}$
Feynman-I-18-4 | $y=\dfrac{m_{1}r_{1}+m_{2}r_{2}}{m_{1}+m_{2}}$
Feynman-I-38-12 | $y=\dfrac{\epsilon h^{2}}{\pi mq^{2}}$
Feynman-I-39-22 | $y=\dfrac{Tkbn}{V}$
Feynman-I-40-1 | $y=n_{0}\operatorname{exp}\left({-\dfrac{gmx}{Tkb}}\right)$
Feynman-I-43-31 | $y=Tkbmob$
Feynman-I-8-14 | $y=\sqrt{\left(-x_{1}+x_{2}\right)^{2}+\left(-y_{1}+y_{2}\right)^{2}}$
Feynman-I-9-18 | $y=\dfrac{Gm_{1}m_{2}}{\left(-x_{1}+x_{2}\right)^{2}+\left(-y_{1}+y_{2}\right)^{2}+\left(-z_{1}+z_{2}\right)^{2}}$
Feynman-test-17 | $y=\dfrac{m^{2}\omega^{2}x^{2}\left(\tfrac{\alpha x}{y}+1\right)+p^{2}}{2m}$
Feynman-test-18 | $y=\dfrac{3\left(H_{G}^{2}+\tfrac{c^{2}k_{f}}{r^{2}}\right)}{8\pi G}$

##  Appendix D Strogatz problems

Table 3 shows the complete Strogatz data set. It can be downloaded from
https://github.com/lacava/ode-strogatz.

Table 3: The Strogatz ODE problem set [La Cava et al., 2016]. Name | Formula
---|---
Bacterial Respiration | $\dot{x}=-\frac{xy}{0.5x^{2}+1}-x+20$
$\dot{y}=-\frac{xy}{0.5x^{2}+1}+10$
Bar Magnets | $\dot{x}=-\sin{\left(x\right)}+0.5\sin{\left(x-y\right)}$
$\dot{y}=-\sin{\left(y\right)}-0.5\sin{\left(x-y\right)}$
Glider | $\dot{x}=-0.05x^{2}-\sin{\left(y\right)}$
$\dot{y}=x-\frac{\cos{\left(y\right)}}{x}$
Lotka-Volterra interspecies dynamics | $\dot{x}=-x^{2}-2xy+3x$
$\dot{y}=-xy-y^{2}+2y$
Predator Prey | $\dot{x}=x\left(-x-\frac{y}{x+1}+4\right)$
$\dot{y}=y\left(\frac{x}{x+1}-0.075y\right)$
Shear Flow | $\dot{x}=\cos{\left(x\right)}\cot{\left(y\right)}$
$\dot{y}=\left(0.1\sin^{2}{\left(y\right)}+\cos^{2}{\left(y\right)}\right)\sin{\left(x\right)}$
van der Pol oscillator | $\dot{x}=-\frac{10x^{3}}{3}+\frac{10x}{3}+10y$
$\dot{y}=-\frac{x}{10}$

##  Appendix E Model parameter search

The success of ParFam depends strongly on a good choice of the model
parameters: The set of base functions $g_{1},...,g_{k}$ and the degrees
$d^{1}_{i}$ and $d^{2}_{i}$, $i\in\\{1,\dots,k+1\\}$, of the numerator and
denominator polynomials of $Q_{1},...,Q_{k+1}$, respectively. On the one hand,
choosing the degrees very small or 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 a global optimization
method. This might prevent ParFam from identifying even very simple functions.

To balance this tradeoff, we allow ParFam to iterate through many different
choices for the hyperparameters describing the model. The user specifies upper
bounds on the degrees $d_{i}^{1}$ and $d_{i}^{2}$ of the polynomials and the
set of base functions $g_{1},\ldots,g_{k}$. 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 procedure is shown in
Algorithm 1. Note that we refer to the rational functions $Q_{1},...,Q_{k}$,
which will be the inputs to the base functions, as the ’input rationals’ and,
therefore, describe the degrees of their polynomials by ’DegInputNumerator’
and ’DegInputDenominator’. Similarly, we denote the degrees of the polynomials
of the output rational function $Q_{k+1}$ by ’DegOutputNumerator’ and
’DegOutputDenominator’.

Input: Maximal Degree Input Numerator $d_{\mathsf{max,in}}^{1}$,

Input: Maximal Degree Output Numerator $d_{\mathsf{max,out}}^{1}$,

Input: Maximal Degree Input Denominator $d_{\mathsf{max,in}}^{2}$,

Input: Maximal Degree Output Denominator $d_{\mathsf{max,out}}^{2}$,

Input: Maximal number of base functions $b_{\mathsf{max}}$

Input: Set of base functions $G_{\mathsf{max}}=\\{g_{1},\dots,g_{k}\\}$.

Output: List of model parameters $\mathcal{L}$ that define the models ParFam
can iterate through.

1 Let $\mathcal{L}=\\{~{}\\}$ be an empty list.

// Start with a polynomial model:

2 $D_{\mathsf{p}}$ = {’DegInputNumerator’: 0, ’DegOutputNumerator’:
$d_{\mathsf{max,out}}^{1}$, ’DegInputDenominator’: 0, ’DegOutputDenominator’:
0, ’baseFunctions’: []}

3$\mathcal{L}$.append($D_{0}$)

// Continue with purely rational models with different degrees:

4 for  _$d^{2}_{\mathsf{out}}=1$ to $d_{\mathsf{max,out}}^{2}$_ do

5 for  _$d^{1}_{\mathsf{out}}=1$ to $d_{\mathsf{max,out}}^{1}$_ do

6 $D_{\mathsf{r}}$ = {’DegInputNumerator’: 0, ’DegOutputNumerator’:
$d^{1}_{\mathsf{out}}$, ’DegInputDenominator’: 0, ’DegOutputDenominator’:
$d^{2}_{\mathsf{out}}$, ’baseFunctions’: []}

7 $\mathcal{L}$.append($D_{\mathsf{r}}$)

8 end for

9

10 end for

11

// Include different combinations of base functions:

12 for  _$b=1$ to $b_{\mathsf{max}}$_ do

13 for  _$d^{2}_{\mathsf{out}}=0$ to $d_{\mathsf{max,out}}^{2}$_ do

14 for  _$d^{1}_{\mathsf{out}}=1$ to $d_{\mathsf{max,out}}^{1}$_ do

15 for  _$d^{2}_{\mathsf{in}}=0$ to $d_{\mathsf{max,in}}^{2}$_ do

16 for  _$d^{1}_{\mathsf{in}}=1$ to $d_{\mathsf{max,in}}^{1}$_ do

17 for  _$B$ as a list with $b$ elements of $G_{\mathsf{max}}$_ do

// Note that base functions can be contained in $B$ multiple times.

18 $D$ = {’DegInputNumerator’: $d^{1}_{\mathsf{in}}$, ’DegOutputNumerator’:
$d^{1}_{\mathsf{out}}$, ’DegInputDenominator’: $d^{2}_{\mathsf{in}}$,
’DegOutputDenominator’: $d^{2}_{\mathsf{out}}$, ’baseFunctions’: $B$}

19 $\mathcal{L}$.append($D$)

20 end for

21

22 end for

23

24 end for

25

26 end for

27

28 end for

29

30 end for

31

return $\mathcal{L}$

Algorithm 1 Traversal of the model parameters

This strategy is comparable to the one proposed by Bartlett et al. [2023],
called “Exhaustive Symbolic 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 set of base functions as we do, but without $\cos$. Our
algorithm, however, only needs to search for the parameters of around 500
functions since it covers many at the same time by employing the global
optimization strategy.

Due to this high complexity, Bartlett et al. [2023] state that they merely
concentrate on one-dimensional problems and, thus, could benchmark their
algorithm only on Feynman-I-6-2a ($y=\exp(\theta^{2}/2)/\sqrt{2pi}$), the only
one-dimensional problem from the Feynman data set [Udrescu and Tegmark, 2020].
This example shows the benefit of employing global search in the parameter
space: 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).

##  Appendix F Hyperparameter settings SRBench ground-truth problems

The hyperparamater settings for the SRBench ground-truth problems are
summarized in Table 4.

Table 4: The model and optimization parameters for the SRBench ground-truth problems Model parameters | Maximal Degree Input Numerator | 2
---|---|---
Maximal Degree Input Denominator | 2
Maximal Degree Output Numerator | 4
Maximal Degree Input Denominator | 3
Base functions |  $\sqrt{}$, $\cos$, $\exp$
Maximal potence of any variable | 3
|  (i.e., $x_{1}^{4}$ is excluded but $x_{1}^{3}x_{2}$ is allowed) |
Optimization parameters | Global optimizer | Basin-hopping
Local optimizer | BFGS
Maximal number of iterations global optimizer | 10
Maximal data set length | 1000
Regularization parameter $\lambda$ | 0.001
Maximal runtime | 8 CPU hours
Maximal number of evaluations | 1,000,000

##  Appendix G Additional plots for the SRBench ground-truth results

Figure 5: Symbolic recovery rate on both SRBench ground-truth data sets
separated. Figure 6: Accuracy solution rate (percentage of data sets with
$R^{2}\\!>\\!0.999$ for the test set) on the SRBench ground-truth problems
separately. Figure 7: Results on the SRBench ground-truth problems. Points
indicate the median test set performance on all problems. The $R^{2}$ Test for
AIFeynman is missing on the plot since SRBench used a higher precision data
type, such that AIFeynman achieved a median $R^{2}$ greater than $1-10^{-16}$.

##  Appendix H Synthetic data set

The synthetic data set for the training and evaluation of DL-ParFam in Section
3.2 is sampled in the following way. We first create a parametric family
$(f_{\theta})_{\theta\in\mathbb{R}^{m}}$ for fixed model hyperparameters,
i.e., for a specific choice of base functions and maximal degrees. For the
experiments, we choose the degree of all denominator polynomials to be 0,
i.e., $Q_{1},...,Q_{k+1}$ are simple polynomials, and their maximal degree
will be set to $2$. As the base functions, we use $\sin$ and $\sqrt{}$.

Now, we aim to sample sparse parameters $\theta\in\mathbb{R}^{m}$ to obtain
interpretable functions. To achieve this, we randomly choose the number of
non-zero coefficients of $Q_{k+1}$ to be $1$ or $2$. Next, we choose those
functions among $Q_{1},...,Q_{k}$ which are used by $Q_{k+1}$ and randomly
select the non-zero coefficients of these. If only one coefficient is chosen,
we ensure that it does not correspond to the constant term. For each non-zero
coefficient chosen that way, we sample a parameter
$\theta_{i}\sim\mathcal{N}(0,9)$.

This way we sample a set of parameters $(\theta^{j})_{j=1,...,K}$ describing
(interpretable) functions $f_{\theta^{j}}$. These functions are evaluated on
the fixed grid $x=-10,-9.9,-9.8,\ldots,9.9$ to create the data sets
$((x_{i})_{i=1,...,m},(y_{i}^{j}=f_{\theta^{j}}(x_{i}))_{i=1,...,m})_{j=1,...,K}$.
Some examples for $f_{\theta^{j}}$ are presented in Table 5.

Table 5: Example formulas from the synthetic data set used to train and
evaluate DL-ParFam in Section 3.2. All coefficients are rounded to three
decimal places.

| Formula
---|---
1. | $y=1.357\sin{\left(4.072x^{2}+1.044\right)}$
2. | $y=-2.909x\sin{\left(8.746x^{2}-1.637x+0.72\right)}+2.591$
3. | $y=4.131x\sin{\left(1.933x^{2}-0.549x+3.205\right)}-2.847\sin{\left(1.933x^{2}-0.549x+3.205\right)}\sqrt{\left|{1.344x+1.678}\right|}$
4. | $y=0.771x\sqrt{\left|{3.293x^{2}+0.878x+1.837}\right|}+3.64\sin{\left(0.824x^{2}-8.78x+1.936\right)}$
5. | $y=-1.375x$
6. | $y=-6.339x\sin{\left(6.961x\right)}+2.891$
7. | $y=0.944x\sin{\left(2.042x^{2}+5.451x+3.97\right)}+0.548x$
8. | $y=-3.907x\sin{\left(5.384x^{2}\right)}+2.681\sin{\left(5.384x^{2}\right)}\sqrt{\left|{0.276x^{2}+2.406x-1.149}\right|}$
9. | $y=4.276x+2.025\sin{\left(3.93x\right)}$
10. | $y=-3.007x\sqrt{\left|{x^{2}}\right|}-2.751\sin{\left(0.231x^{2}-1.3\right)}$
11. | $y=-2.189x^{2}-2.828\sin{\left(0.188x^{2}+0.63\right)}$
12. | $y=-0.317\sin{\left(-2.351x^{2}+1.448x+2.344\right)}-6.667$
13. | $y=2.154x+3.064\sin{\left(4.773x^{2}\right)}\sqrt{\left|{4.491x-0.423}\right|}$
14. | $y=-0.772x^{2}+0.333x\sin{\left(2.938x^{2}+2.245x\right)}$
15. | $y=-6.929\sin{\left(2.451x^{2}+2.37x+4.415\right)}\sqrt{\left|{2.486x-1.428}\right|}-1.873\sin{\left(2.451x^{2}+2.37x+4.415\right)}$
16. | $y=-0.528x^{2}-0.601\sin{\left(3.249x^{2}\right)}$
17. | $y=-3.77\sin{\left(0.577x\right)}\sqrt{\left|{0.981x-2.192}\right|}$
18. | $y=1.076\sin{\left(0.189x\right)}$
19. | $y=-2.046x$
20. | $y=1.608\sin{\left(5.982x-2.644\right)}\sqrt{\left|{x^{2}}\right|}$

##  Appendix I Hyperparameter settings synthetic benchmark

The hyperparamater settings for the experiments on the synthetic benchmark in
Section 3.2 are summarized in Table 6.

Table 6: The model and optimization parameters for the benchmark with the synthetic data sets. Data set parameters | Training set size | 2,000,000
---|---|---
Validation set size | 10,000
Test set size (for the NN) | 10,000
Test set size (for DL-ParFam vs ParFam) | 100
Maximal number of non-zero coefficients of $Q_{k+1}$ | 2
Model parameters (data creation and training) | Maximal Degree Input Numerator | 2
Degree Input Denominator | 0
Degree Output Numerator | 2
Degree Input Denominator | 0
Base functions |  $\sqrt{}$, $\sin$
Maximal potence of any variable | 2
Optimization parameters (Neural network pre-training) | Optimizer | ADAM
Loss | BCE
Number epochs | 20,000
Number batches | 500
Number hidden layers | 1
Number hidden neurons per layer | 200
Learning rate | 0.0001
Optimization parameters (ParFam and DL-ParFam) | Global optimizer | Basin-hopping
Local optimizer | BFGS
Maximal number of iterations global optimizer | 10, 20, 50, 100
Regularization parameter $\lambda$ | 0.001
Maximal runtime | no limit
Maximal number of evaluations | no limit