GSR: A Generalized Symbolic Regression Approach
Authors: Tony Tohme
Dehong Liu
KAMAL YOUCEF-TOUMI
Venue: Accepted by TMLR ()

Abstract: Identifying the mathematical relationships that best describe a dataset remains a very challenging problem in machine learning, and is known as Symbolic Regression (SR). In contrast to neural networks which are often treated as black boxes, SR attempts to gain insight into the underlying relationships between the independent variables and the target variable of a given dataset by assembling analytical functions. In this paper, we present GSR, a Generalized Symbolic Regression approach, by modifying the conventional SR optimization problem formulation, while keeping the main SR objective intact. In GSR, we infer mathematical relationships between the independent variables and some transformation of the target variable. We constrain our search space to a weighted sum of basis functions, and propose a genetic programming approach with a matrix-based encoding scheme. We show that our GSR method is competitive with strong SR benchmark methods, achieving promising experimental performance on the well-known SR benchmark problem sets. Finally, we highlight the strengths of GSR by introducing SymSet, a new SR benchmark set which is more challenging relative to the existing benchmarks.

Published in Transactions on Machine Learning Research (12/2022)
GSR: A Generalized Symbolic Regression Approach
Tony Tohme
tohme edu
Massachusetts Institute of Technology
Dehong Liu
Gudh@merr..c com
Mitsubishi Electric Research Laboratories
Kamal Youcef- Toumi
youcef@mil. edu
Massachusetts Institute of Technology
// Reviewed on OpenReview: https: // optpr/veew.eee.... net/ forum? id= tDuU
Abstract
Identifying the mathematical relationships that best describe dataset remains a very chal- lenging problem in machine learning, and is known as Symbolic Regression (SR). In contrast to neural networks which are often treated as black boxes, SR attempts to gain insight into the underlying relationships between the independent variables and the target variable a given dataset by assembling analytical functions. In this paper, we present GSR, a General- ized Symbolic Regression approach, by modifying the conventional SR optimization problem formulation, while keeping the main SR objective intact. In GSR, we infer mathematical relationships between the independent variables and some transformation of the target vari- able. We constrain our search space to weighted sum of basis functions, and propose genetic programming approach with matrix-based encoding scheme. We show that our GSR method is competitive with strong SR benchmark methods, achieving promising experimental performance on the well-known SR benchmark problem sets. Finally, we high- light the strengths of GSR by introducing SymSet, new SR benchmark set which more challenging relative to the existing benchmarks.
1 Introduction
Symbolic regression (SR) aims to find mathematical expression that best describes the relationship between the independent variables and the largel. (or dependent. variable based on given dataset. By inspecting the resulting expression, may be able to identify nontrivial relations and/or physical laws which can provide more insight into the system represented by the given datasct. SR has gained tremendous interest and attention from researchers over the years for many reasons. First, many rules and laws in natural sciences (e.g. in physical and dynamical systems (Schmidt & Lipson 2009 Quade et al. 2016)) are accurately represented by simple analytical equations (which can be explicit (Brunton et al. 2016)) or implicit (Mangan et al.) 2016 Kaheman et al.) 2020)). Second, in contrast to neural networks that involve complex input- output mapping, and hence often treated black boxes which are difficult to interpret, SR is very concise and interpretable. Finally, symbolic equations may outperform neural networks in out-o--distributio gencralization (especially for physical problems) (Cranmer et al. 2020).
SR does not require priori specification of a model. Conventional regression methods such as least squares
(Wild & Seber 1989), likelihood- baaeed (Edwards) 1984 Pawitan 2001 Tohme et al. 2021), and Bayesian
regression techniques [Lee) 1997 Leonard Hsu 2001 Tohme 2020 Tohme et al., 2020 Vanslette et al.,
2020) use fixed-form parametric models and optimize for the model parameters only. SR seeks find both model structure and its associated parameters simultaneously.

---

Published in Transactions on Machine Learning Research (12/2022)
Related Work. The SR problem has been widely studied in the literature Orzechowski al.
2018
La Cava
et al. 2021 ). SR can be a very challenging problem and is thought to be NP-hard (Lu et al. 2016 Udrescu
& Tegmarkk 2020 Petersen et al.) 2021; Virgolin & Pissis 2022]. It can also be computationally expensive as the scarch space is very wide (or complex) containing expressions of any size and length de França (2018), this issue being exacerbated with the dimension of the input feature vector (i.e. the number of
independent variables). Several approaches have been suggested over the years. Most of the methods
use genetic (or evolutionary) algorithms (Koza & Koza 1992 Schmidt & Lipson
2009
Bäck al., 2018;
Virgolin et al. 2019). Some more recent methods are Bayesian in nature (Jin al.
2019
some are physics-
inspired (Udrescu & Tegmark 2020), and others use divide-and-conquee (Luo et al.
2017
and block building
algorithrms (Chen et a 7b1 2018a b). Lately, researchers proposed using machine learning algorithms
and neural networks to solve the SR problem [Martius & Lampert
2016
Sahoo et al.) 2018
Udrescu el al.
2020 Ahn et al. 2020 Al-Roomi & El-Hawary 2020
Kim et al. 2020
Kommenda et al.
2020
Burlacu
ct al. 2020 Biggio ct al. 2021 Mundhenk ct al., 2021
Petersen et al.,
2021
Valipour ct al.,
2021
Razavi
& Gamazon 2022 Zhang el al. 2022a] d Ascooi et al.
2022 Kamienny
al.
2022; Zhang al.
2022b).
Furthermore, some works suggested constraining the search space of functions to generalized linear space
Nelder & Wedderburn 1972) (e.g. Fast Function eXtraction (McConaghy 2011), Elite Bases Regression
Chen et al.) 2017ap, etc.) which proved to accelerate the convergence of genetic algorithms significantly (at the expense of sometimes losing the generality of the solution (Tuo et al. (017)).
Most of the SR methods use a tree- based implementation, where analytical functions are represented (or encoded) by expression trees. Some approaches suggested encoding functions as an integer string (ONNiil & Ryan 2001), others proposed representing them using matrices (Tuo & Zhang 2012 Chen et al., 2017a
de França, 2018 de Franca & Aldeia 2020). As we will discuss in later scctions, implementation relies matrices to encode expressions.
Our Contribution. We present Generalized Symbolic Regression (GSR), by modifying the conventional SR optimization problem formulation, while keeping the main SR objective intact. In GSR, we identify mathematical relationships between the independent variables (or features) and some transformation of the target variable. In other words, learn the mapping from the feature space to transformed target space (where the transformation applied to the target variable is also learned during this process). To find the appropriate functions (or transformations) to be applied to the features as well as to the targets, we constrain our search space to weighted sum of basis functions. In contrast to conventional tree- based genetic programming approaches, we propose a matrix-based encoding schome to represent the basis functions (and hence the full mathematical expressions). We run a series of numerical experiments on the well-known SR benchmark datasets and show that our proposed method is competitive with many strong SR methods. Finally, we introduce SymSet, new SR benchmark problem set that is more challenging than existing benchmarks.
2 Notation and Problem Formulation
Consider the following regression task. We are given a dataset D= = {x,, Ji consisting of N i..dd paired examples, where x; TRd denotes the 7th d-dimensional input feature vector and yi R represents the corresponding continuous target variable. The goal of SR is to search the space of all possible mathematical expressions S defined by set of given mathematical functions (e.g., exp, In, sin, cos) and arithmetic operations (e.g., +, -, x, +), along with the following optimization problem:
arg min f
Σ
[f(x;) ya
(1)
JES
where f is the model function and f is the optimal model.
2

---

Published in Transactions on Machine Learning Research (12/2022)
3
Generalized Symbolic Regression (GSR)
In this section, we introduce our Generalized Symbolic Regression (GSR) approach. We present its problem formulation, and discuss its solution and implementation.
3.1 Modifying the goal of symbolic regression
As highlighted in Section 2 the goal of SR is search the function space to find the model that best fits the mapping between the independent variables and the target variable (i.e. the mapping between X; and yt, for all i). Since the main objective of SR is to recognize correlations and find non-trivial interpretable models (rather than making direct predictions), we modify the goal of SR; we instead search the function space to find the model that best describes the mapping between the independent variables and a transformation of the target variable (i.e. the mapping between x; and some transformation or function of yi, for all i), Formally, propose modifying the goal of SR to scarch for appropriate (model) functions from space of all possible mathematical expressions S defined by set of given mathematical functions (e.g., exp, In, sin, cos) and arithmetic operations (c.g., +, - x, +), which can be described by the following optimization problem:
f',g" arg min fur ES
E [f(x,) (((x)))9(ux) g(v))2
where f and g" are the optimal analytical functions. In other words, instead of searching for mathematical expressions of the form y = fxx as is usually done in the SR literature, the proposed GSR approach attempts to find expressions of the form (vy) = fxx). We illustrate this concept in Table
Table 1: GSR finds analytical expressions of the form g(v) fxx instead of y f(x).
Ground Truth Expression
Learned Expression
y Vx +
y² x +
y 1/(3x) + +x))
y 1 3x1 + x
y = (2xr + +x22
In(y) In + x2)
У In(x} + 42)22)
el =² -4x102
= " cos3))
In(y) = x9 + 2x² + cos(ass
Although the main goal of GSR is to find expressions of the form g f(x), we may encounter situations where it is best to simply learn expressions of the form y f(x) (i.c. g(y) y)- For instance, consider the ground truth expression
sin(æ1) + 2zz. In this case, expect to learn the expression exactly as is (i.e. (()) = y and f(x) = sin(r1) + 2x2) as long as the right basis functions (i.e. sin(x) and x in this case) are within the search space, as we will see in the next sections.
Making predictions. Given a new input feature vector x., predicting y, with GSR is simply matter of solving the equation s(y) fxx) for y, or equivalently, g(y) f(x.) 0. Note that f(x.) is known quantity and is the only unknown. If g(-) is invertible function, then ys be casily found using y. ((xx)) If g(-) is not invertible, then y, will be the root of the function h(y) g(v) - f(x.). Root-finding algorithms include Newton's method. Whether the function 9(.) is invertible or not, we might end up with many solutions for ys (an invertible function, which is not one-to-one, can lead to more than one solution). In this case, we choose y, be the solution that belongs to the range of y which can be determined from the training dataset.
3

---

Published in Transactions on Machine Learning Research (12/2022)
3.2 A new problem formulation for symbolic regression
Now that we have presented the goal of proposed GSR approach (summarized by Equation 2). need to constrain the search space of functions S to reduce the computational challenges and accelerate the convergence of our algorit Inspired by McConaghy (2011 : Chen et al. (2017a) as well as classical system identification methods (Brunton et al. 2016), we confine S to generalized linear models, i.c. to functions that can be expressed as a linear combination (or as a weighted sum) of basis functions (which can linear or nonlinear). In mathematical terms, for given input feature vector x, and corresponding target variable Vi, the search space is constrained to model functions of the form:
Mr
Me
f(x)) - Σ (x)) 1=1
glyb - Σ (,0,(1)) ==1
(3)
where d)(+) and (-y(-) are the basis functions applied to the feature vector x; and the target variable Vi, respectively, Mg and My denote the corresponding mumber of basis functions involved, respectively. In matrix form, the minimization problem described in Equation 2 is equivalent to finding the vectors of
coefficients a = [ai ............. "aa.]" and B = [3, ... SM,,] such that:
a B = arg min |Xx - 1 a,
(4)
where
(x1)
(x2(x)) ... SN, (x1)
(i((1) (2((1)
VMs, (v1)
X
Y
(5)
(xxN
(x())
(xN)
s (NN) (2() UM, (yN)
Note that we examine the minimization problem as expressed in Equation we can indeed minimize |Xa Y by simply setting ax and 8 0 which will not lead to meaningful solution to GSR problem. In addition, to avoid reaching overly complex mathematical expressions for f(-) and g(-), we are interested in finding sparse solutions for the weight vectors a and 3 consisting mainly of zeros which results in simple analytical functions containing only the surviving basis functions (i.e. whose corresponding weights are nonzero). This is closely related to sparse identification of nonlinear dynamics (SINDy) methods [Buntoon et al., 2016). To this end, we apply L₁ regularization, also known as Lasso regression Tibshirani 1996), by adding a penalty on the L1 norm of the weights vector (i.e. the sum of its absolute values) which leads to sparse solutions with few nonzero coefficients. In terms of our GSR method, Lasso regression automatically performs basis functions selection from the set of basis functions that under consideration.
Putting the pieces together, we reformulate the minimization problem in Equation as a constrained Lasso regression optimization problem defined as
w" = arg min ||A 2 + A||W||
(6)
s.t. || |wll2 =
where À > 0 is the regularization parameter, and
a
x
YJ.
w
[a
ams B1 B.....
(7)
B
3.3 Solving the GSR problem
To solve the GSR problem, first present our approach for solving the constrained Lasso problem in Equation assuming some particular sets of basis functions given. We then outline our genetic program- ming (GP) procedure for finding the appropriate (or optimal) sets of these basis functions, before discussing our matrix-based based encoding scheme (to represent the basis functions) that we will use in our GP algorithm.
4

---

Published in Transactions on Machine Learning Research (12/2022)
Algorithm 1: Solving the constrained Lasso optimization problem using ADMM
Input: A, A: P. Wo: 20, ug
Output:
function SOLVEADMM(A. X, Ps wu, 2u. Ug Initialization: w + wp, 20: u M40: while Not Converge do
w - (2ATA + pI) p (x u): w/||wl12:
z - Sxip w u); u - u + w z;
end end function
3.3.1 Solving the Lasso optimization problem given particular sets of basis functions
We assume for now that, in addition to the datasct D (x, v we are also given the sets of basis functions {x)(x)) Me and (1)) Afe used with the imput feature vector and its corresponding target variable Vi, respectively, for all i. other words, we assume for now that the matrix A in Equation 7 is formed based on particular sets of basis functions (i.e. (x)(x:)) Ms
j=1 and (w))), and we are mainly interested in solving the constrained optimization problem in Equation o Applying the alternating direction method of multiplicrs (ADMM) [Boyd et al. 2011). ) the optimization problem in Equation 6 can be written as
w" = arg min ||w + A| 1|1
s.t. |
(8)
w
0
where À > 0 is the regularization parameter. The scaled form of ADMM (see Boyd et al. 2011 for details) for this problem is
Wy < arg min C,(W, Zk-1, u6-1) |=|1==1
(9)
2 < S2/o( (wk ux=11)
<- u
where u the scaled dual vector, and
Cow x, u) = |Aw + w - 2 + =u12
(10)
where p > 0 is the penalty parameter and the soft thresholding operator defined as
a -- K a > K
Sx(a) = 0 la| <
(11)
a+k
a
< -K
To find the minimizer we the first step of the ADMM algorithm above (in Equation 9 we first com- pute the gradient of the function Cp(w, (w, 2k- ux-1) with respect to w, set it to zero, and then normalize the resulting vector solution:
WnL,(W, zk (6-11
(12)
-2ATAA Aun + + (wk 2k-1 + 2k-1

---

Published in Transactions on Machine Learning Research (12/2022)
It follows that
= (2.* A pl) p1) -p ux-1)
(13)
Wk = ww w
Aggorthm@ooutlnge the overall process for solving the constrained Lasso optimization problem in Equation b for given matrix A, regularization parameter X, penalty parameter μ, and initial guesses Wo, 20, u0-
3.3.2 Finding the appropriate sets of basis functions using genetic programming
Now that we have presented Algorithm that solves the constrained Lasso optimization problem Equation 0 for particular sets of basis functions, we go through our procedure for finding the optimal sets of basis functions (and hence, the optimal analytical functions f() and 9(-)).
Encoding Scheme
Most of the SR methods rely on expression trees in their implementation. That is, each mathematical expression represented by a tree where nodes (including the root) encode arithmetic operations (e.g. + X, +) or mathematical functions (e.g. cos, sin, exp, In). and leaves contain the independent variables (i.c. 21...., Faa) or constants. Inspired by Luo & Zhang (2012;; Chen et al. (2017a),, use matrices instead of trees to represent the basis functions. However, propose general enough to handle/recover wide range of expressions.
encoding scheme that we believe is
We introduce the basis matrices B and B to represent the basis functions ((-) and (b(+) used with the feature vector x and the target variable y, respectively. The basis matrices 8 and B are sizes ng X Xm and gg 1 respectively (i.c. B is a column vector), and take the form
b m g
1,1
B
B
(14)
6
b
by
12 BB..
where the entries b and 1 are all integers. The first column b", of Bdd and the first (and only) column b of B indicate the mathematical function (or transformation) to be applied (on the input feature vector x and the target variable y, respectively). In Bd, the second column 6",2 specifies the type of argument (see Table 2), and the remaining n. m= -2 columns be "........ 6%. me indicate which independent variables (or features) are involved (i.e. the active operands). The quantity the represents the maximum total multiplicity of all the independent variables included in the argument. Note that ng and ng specify the number of transformations Lo be multiplied logether (i.e. each basis function o(-) and (>)- will be a product of n and B transformations, respectively).
table of mapping rules that is very straightforward understand and employ. For instance, consider the mapping rules outlined in Table 2 where d is the dimen- sion of the input feature vector. As we will see in our numerical experiments in Section 4 we will adopt this table for many SR benchmark problems. Other mapping tables are defined according to different benchmark problems (more details about the SR benchmark problem specifications can be found in Appendix C The encoding from the analytical form of a basis function to the basis matrix is straightforward. For example, for d = 3, nB = and mB = 5 (i.e. n, = 3), the basis function Φ(x) = æ2 cos(afa2) ln(2 + xx) can be generated according the encoding steps shown in Table 3
The mncodin//ecooding /ecodiing process happens according
Each SR. benchmark problem specific (or library) of allowable mathematical functions (e.g. cos, sin, exp, log): and hence, we mainly modify the first rows of the mapping tahles.

---

Published in Transactions on Machine Learning Research (12/2022)
Example table of mapping rules for basis function. The identity operator is denoted by o!.
Table 2:
Table 3: Encoding steps corresponding to the
basis function ((x) 22 c2) In(x1 + xz).
b
Step
T
arg
01
92
1/
Update
Transformation (2)
e
co8
sin
exp
a
x2
Ti Ti(x) 22
b,,2
2
cos
3:2
T2(x) = cos (xf
Argument Type (arg)
Σ II
In
23 skip
Ta(x) = In(x1 23)
Ta(x)
be.3.'.. bamm
1
...
Variable (v)
skip
Final Update: 6(x) (x)) (x)x) Tx(x) (x))
Based on the mapping rules in Table 2 and the encoding steps in Table o(x) 22 cos()f (n(21 + 2z) can be described by x 5 matrix as follows:
3 the basis function
1
0
2
B
2
2
1
1
2
(15)
5
1
1
3
0
Remark 3.1. In Table denotes entries that are ignored during the construction of the basis function. The argument type in Step is I which implies that we only select the first variable (encoded by b.,3) out of the n, variables as an argument, and hence the entries corresponding to 82 and 03 are ignored. Similarly, the transformation in Step is == which implics that the argument type and the n. variables are all ignored. These the only cases where some entries are ignored during the construction process. The same ignored entries are reflected in the matrix 8 using More encoding examples can be found in Appendix B
Remark 3.2. The term 'skip' can be thought of as 0 or 1 when the argument type is summation Σ multiplication II respectively. To account for the case where the argument type is x, we let b,,3 E {l,...,d d (i.c. we exclude 0) as b,, is the only entry considered in this case (see Remark 3.1
Remark 3.3. The same basis function can be represented by several matrices for three reasons: i) Each basis function product of transformations where each transformation is represented by row in the basis matrix. Hence, a basis matrix for the same basis function is formed by simply swapping rows. ii) When the argument type is Σ or II, the order of the variables (including 'skip') starting from the third column of the matrix 8 does not affect the expression. Hence a new basis matrix for the same basis function formed by simply swapping these columns.
iil) As mentioned in Remark 3.1 some entrics are ignored in some cases. Hence a new basis matrix for the same basis function formed by simply modifying these entries.
Remark: 3.4. In the example above, we showed how we can produce the matrix B to represent basis func- tion @(x). A similar (and even simpler) procedure be applied to produce the matrix 8 that represents basis function (w);; we only need mapping table corresponding to the set of allowable transformations (e.g. the first two rows of Table 2).
Note that the decoding from the basis matrix to the expression of basis function is trivial; we go through the rows of the basis matrix and convert them into transformations according to mapping ta- ble (e.g. Table 2). before finally multiplying them together. Also note that the search space of basis functions (mainly ((-)) huge in general which makes enumeration impractical, and hence, we will rely GP for effective search process.

---

Published in Transactions on Machine Learning Research (12/2022)
Genetic Programming (Evolutionary Algorithm)
The SR problem has been extensively studied in the literature, and a wide variety of methods has been suggested over the years to tackle it. Most of these methods are based on genetic programming (GP) (Koza
& Koza
1992 Schmidt & Lipson 2009 Bäck et al.) 2018) Virgolin et al. 2019). This is a houristic search technique that tries to find the optimal mathematical expression (in the SR context) among all possible expressions within the search space. The optimal (or best) expression is found by minimizing some objective
function, known as the fitness function.
GP is an evolutionary algorithm that solves the SR. problem. It starts with an initial population (or first generation) of Np randomly generated individuals (i.c. mathematical expressions), then recursively applics the selection, reproduction (crossover), and mulation operations until termination. During the selection operation, the fitness of each of the Np individuals of the current generation is evaluated (according to the fitness function), and the ny fittest (or best) individuals are selected for reproduction and mutation (the selected individuals are part of the new generation and can be thought of as parents). The reproduc- tion (crossover) operation generates new individuals (offsprings) by combining random parts of two parent individuals. The mutation operation produces new individual by changing random part of some par- ent individual. Finally, the recursion terminates, when some individual reaches a predefined fitness level (i.e. until some stopping criterion is satisfied).
In our GSR approach, we use a slightly modified version of the GP algorithm described above. Each individual in the population initially consists of two sets of Ms and Ms randomly generated basis functions encoded by basis matrices. Such matrices will form the functions f(.) and g(-) to be used with the input feature vector x and the target variable y, respectively. This is different from the GP algorithm described above where individuals typically represent the full mathematical expression or function as a whole. In addition, the N₂ individuals in the population of a new generation consist of the n, fittest individuals of the current generation in addition to N2 - np individuals generated as follows. With probability 1. a new individual is generated (reproduced) by randomly combining basis functions (i.c. basis matricos) from two parent individuals (i.e. crossover) selected from the "p surviving individuals. With probability 4, new individual is generated by randomly choosing one of the ny surviving individuals, and replacing (mutating) some of its basis functions (i.e. basis matrices) with completely new ones (i.e. randomly generated). With probability 2: a completely new individual is randomly generated (in the same way we generate the individuals of the initial population). Randomly generating individuals enhances diversity in the basis functions and avoids reaching a plateau. Indeed, this is just one of many ways that can be followed to apply some sort of crossover crosoovrr/muaatin /mutation on individuals defined by their sets of basis functions instead of their full mathematical expression. A pseudocode of our proposed GSR algorithm is provided in Appendix A
I
4 Experimental Results
We evaluate our proposed GSR method through a series of numerical experiments on a number of common SR benchmark datasets. In particular, we compare our approach to existing state-of-the-art methods using three popular SR benchmark problem sets: Nguyen (Uy et al. 2011). .Jin (Jin et al. 2019), and Neat (Truiilo) let al., 20166). In addition, we demonstrate the benefits of our proposed method on the recently introducod SR benchmark dataset called Livermore (Mundhenk et al., 2021), , which covers problems with a wider range of difficulty compared to the other benchmarks. Finally, we introduce new and more challenging set of SR. benchmark problems, which we call SymSet, mainly for two reasons: i) Our GSR algorithm achieves perfect scores on Nguyen, and almost perfect scores on Jin, Neat, and Livermore, and hence we introduced a benchmark problem set that is more challenging, ii) The existing SR benchmark problem sets do not really reflect the strengths of our proposed method, and thus we designed SymSet to explicitly highlight the benefits we gain from using our proposod approach. SymSet contains benchmark problems with similar properties as Nguyen, Jin, Neat, and Livermore benchmarks, but with an additional function composition (or symbolic layer). Fach SR benchmark problem consists of ground truth expression, a training and test dataset, and set (or libary) of allowable arithmetic operations and mathematical functions. Specifications of all the SR benchmark problems are described in Appendix Hyperparameters and additional experiment details are provided in Appendix A
8

---

Published in Transactions on Machine Learning Research (12/2022)
Table 4: Recovery rate comparison of GSR against several algorithms on the Nguyen benchmark set over 100 independent runs. The formulas for those benchmarks are shown in Appendix Table 19
Recovery Rate (%)
Benchmark
Expression
GSR
NGGPPS
DSR
Eureqa
Nguyen-1
x
xx²
+
+ I
100
100
100
100
Nguyen-2
y +
x
+ +x² + tx +
100
100
100
100
Nguyen-3
= xº + x
+ +
100
100
100
95
Nguyen-4
γ x +
x
+
+ x + x
100
100
100
70
Nguyen-5
sin(x²)
cos (x) -
100
100
72
73
N
y sin(x) sin(2 Ø x²))
100
100
100
100
Nguyen-7
= In(2 + 1) + In +1) Ч VI
100
97
35
85
Nguyen-8
100
100
96
Nguyen-9
sin(x1) ) + sin 2sin(z1) cos(
100
100
100
100
Nguyen-10
100
100
100
64
Nguyen-11
y 2
100
100
100
100
Average
100
99.73
91.18
80.64
Across our experiments, we compare our GSR approach against several strong SR benchmark methods: Neural-guided genetic programming population sceding (NGGPPS): A hybrid approach of neural- guided search and GP, which uses a recurrent neural network (RNN) to seed the starting population for GP (Mundhenk et al., p.21p. NGGPPS achieves strong results on the well-known SR benchmarks.
Deep Symbolic Regression (DSR): A reinforcement learning method that proposes a risk-seeking policy gradient to train RNN to produce better-fitting expressions (Petersen et al., 2021). DSR is the "RNN only" version of NGGPPS, and is also considered a strong performer on the common SR benchmarks. Bayesian Symbolic Regression (BSR): A Bayesian framework which carcfully designs prior dis- tributions to incorporate domain knowledge (e.g. preference of basis functions or tree structure), and which employs efficient Markov Chain Monte Carlo (MCMC) methods to sample symbolic trees from the postcrior distributions qJin et al.] 2019).
Neat-GP: GP approach which uses the NeuroEEvlution of Augmenting Topologies (NEAT) algorithm that greatly reduces the effects of bloat (i.e. controls the growth in program size) (Trujillo et al. 2016). PSTree: A piece-wise non-lineat SR method based decision tree and GP techniques (Zhang et al. 2022a
PSTree can generate explainable models with high accuracy in short period of time. PSTree is the cur- rent top performer on SRBench datasets (La Cava et al. 2021), achieving stale-of-the-arl performance and
beating other competitive SR methods such as Operon Kommenda et al. 2020 Burlacu et al. 2020 and AI Feynman [Udrescu et al.) 2222)).
PySR: A fast and parallelized SR method in Python/Julia (Cranmer 2020 which uses evolutionary algo-
rithms to search for symbolic expressions by optimizing a particular objective; the metric used for scoring equations is based on the work by Cranmer et al. (2020).
gplearn: A Koza-style SR method in Python, which starts with a random population of models, and then iteratively performs tournament selection, crossover, and mutation (Koza & Koza 1992).
We first compare GSR against NGGPPS, DSR, well as Eureqa (a popular GP-based commercial software proposed in Schmidt & Lipson (2009) the Nguyen benchmarks. We follow their experimental procedure and report the results in Table We use recovery rate as our performance metric, defined as the fraction of independent training runs in which an algorithm's 's resulting expression achieves exact symbolic equivalence compared to the ground truth expression (as verified using computer algebra system such as SymPy (Meurer) et al. 2017)). Table 4 shows that GSR significantly outperforms DSR and Eureqa in exactly recovering the Nguyen benchmark expressions. As NGGPPS achieves nearly perfect scores on the Nguyen benchmarks, GSR shows only slight improvement (on Nguyen-7) compared to NGGPPS. However, GSR exhibits faster runtime than NGGPPS; by running cach benchmark problem, GSR takes average of 2.5 minutes per run on the Nguyen benchmarks compared to 3.2 minutes for NGGPPS. Runtimes on individual Nguyen benchmark problems are shown in Appendix Table II
9

---

Published in Transactions on Machine Learning Research (12/2022)
Table 5: Comparison of mean root-mean- square error (RMSE) for GSR against several methods problem set over 50 independent runs. The formulas for these benchmarks are shown in
Table 6: Comparison of median RMSE for GSR against several methods on the Neat benchmark problem set over 30 independent runs. The formulas for these benchmarks are shown in Appendix Table 19
the Jin benchmark
Appendix Table 19
Median RMSE
Benchmark
GSR
NGGPPS
DSR Neat-GP
Mean RMSE
Neat-1
0.0779
Benchmark
GSR NGGPPS
DSR
BSR
Neat-2
0
0.0576
Neat-3
0
0
0.0041
0.0065
Jin-1
0
0
0.46
2.04
Neat-4
0
0.0189
0.0253
Jin-2
6.84
Neat-5
0.0023
Jin-3
0
0
0.00052
0.21
Neat-6
2.0 10-4
6,1 10-6
0.2378
0.2855
Jin-4
0
0.00014
0.16
Neat-7
0.0521
1.0028
1.0606
1.0541
Jin-5
0
0.66
Neat-8
4.0 x 10-4
0.0228
0.1076
0.1
Jin-6
0.018
0
2.23
4.63
Neat-9
8.1 10-9
0.1511
0.1202
Average
0.0030
0
0.45
2.42
Average
0.0059
0.1139
0.1756
0.1977
We next evaluate GSR on the Jin and Neat benchmark sets. The results are reported in Tables 5 and i respectively. A RMSE value of indicates exact symbolie equivalence. From Table b we can clearly observe that GSR out outperforms DSR and BSR and performs nearly as good as NGGPPS recovering all the Jin problems (accross all independent runs) except Jin-6. Table 6 shows that GSR outperforms all other methods (NGGPPS, DSR, and Ncat-GP) on the Neat benchmarks. Note that expressions containing divisions (i.e. Neal-6, Neat-8, and Neat-9) are not exactly recovered by GSR (i.e. only approximations are recovered) since the division operator is not included in our scheme (see Appendix for details).
We then run experiments on the Livermore benchmark set which contains problems with large range of difficulty. In addition to NGGPPS and DSR, we compare against NGGPPS using the soft length prior (SLP) and hicrarchical entropy regularizer (HER) recently introduced in Larma et al. (2021). We also compare against a recently proposed method, known by genetic expert-guidee learning (GEGL) (Ahn el al. 2020 which trains molecule- generrtingn deep neural network (DNN) guided with genetic exploration. Table shows that our GSR method outperforms all other methods both the Nguyen and Livermore benchmark sets, beating NGGPPS+SLP /HER which was the top performer these two benchmark sets.
We highlight the strengths of GSR on the new SymSet benchmark problem set, and show the benefits of searching for expressions of the form g(v) = f(x) instead of y f(x). Typical expressions, with exact symbolic equivalence, recovered by GSR are shown in Appendix Table 26 The key feature of GSR lies in its ability to recover expressions of the form g = fx). To better highlight the benefits offered by this feature, we disable it by constraining the search space in GSR to expressions of the form y = f(x) (which is the most critical ablation). We refer to this special version of GSR as s- S-G.. Note that most of the SymSet expressions cannot be exactly recovered by s-GSR (i.e. they can only be approximated)). We compare the performance of GSR against s-G s-GSR on the SymSet benchmarks in terms accuracy and runtime (sce Table X The results clearly show that GSR is faster than s-GR,, averaging around 2 minutes per run on the SymSet benchmarks compared to 2.27 minutes for s-GSR (i.e. 11% runtime improvement). In addition, GSR is more accurate than s-GRR by two orders of magnitude. This is due to the fact that GSR exactly recovers the SymSet expressions across most of the runs, while s-GSR only recovers approxinations for most of these expressions. This reflects the superiority of GSR over s-SR, which demonstrates the benefits of learning expressions of the form g(y) = fx) in SR tasks. We further compare GSR against several strong SR methods with similar (or better) expression ability. In particular, we experiment on SymSet with NGGPPS, PSTree, PySR, and gplearn (see Table 3 GSR is more accurate than all these methods by three orders of magnitude, which further demonstrates the advantage of our proposed approach. As for the runtime, PSTree is the fastest method, averaging around 16 seconds per run on the SymSet expressions, while maintaining solid accuracies. This comes as no surprise given its state-of- oftte-of-the-art performance on SRBench datasets (La Cava et al. 2021
10

---

Published in Transactions on Machine Learning Research (12/2022)
Table 7: Recovery rate comparison of GSR against several algorithms on the Nguyen and Livermore benchmark sets over 25 independent runs. Recovery rates on individual benchmark problems are shown in Appendix Table IO
Table 8: Average performance in mean RMSE and runtime, along with their standard errors, for GSR against s-GSR and several strong SR methods on the SymSet benchmark problem sets over 25 independent runs. Mean RMSE and runtime values on individual benchmark problems are shown in Appendix Table 12
Recovery Rate (%)
SymSet. Average Mean RMSE
All
Nguyen Livermore
Runtime (sec)
GSR
90.59
100.00
85.45
GSR
2.66 10 t 1.59
10
120.84
t 4.22
NGGPPS+SLP/HER
82.59
92.00
77.45
GSR 2.56 x 10 1 5.27
10 3
136.19
1 4.56
NGGPPS
78.59
92.33
71.09
PSTree
4.57 x 10
1 7.98
10-2
16.23
0.67
NGGPPS
4.65 x 10
1 1.244
10-1
158.57
2.59
GEGL
66,
86.00
56.36
PySR
4.99 x 10
: 1.76
10-1
87.07
21.6
DSR
49.18
83.58
30.41
gplearn
7.22 x 10
1.64 10-1
163.86
2.94
5 Discussion
Limitations. GSR, including state-of-the-art methods, have difficulty with expressions containing divisions. For GSR, this is due to the way we define our encoding scheme. Other methods fail even though the division is included in their framework. GSR can overcome this issue by modifying its encoding scheme to include divisions within the basis functions (at the expense of significantly increasing the complexity of the search space). Another limiting factor to GSR is that it cannot recover expressions containing composition of functions, such as y ex0((2) + In(x). This could be overcome by modifying the search space (e.g. one could expand the definition of basis function to account for composition of functions up to some number of layers, or completely modify the search space to symbolic neural network as in Martius & Lampert (2016); Sahoo et al. (2018); Kim et al. (2020)) Another challenging task for GSR is to reach, although expressible, expressions containing multiple complex basis functions simultaneously. This can be due to the choice of the hyperparameters or the GP search process. A more elaborate discussion about the limitations of GSR can found in Appendix E These limitations will addressed in future paper. Indeed, there are plenty of expressions that still cannot be fully recovered by GSR. This is the case for all other SR methods as well.
Closely related work. There has been growing attention the SR task with non explicit (or implicit) mathematical equations and several works have been attempted to address this interesting task. In partic- ular, implicit sparse identification of nonlinear dynamics (implicit-SINDy) (Mangan et al. 2016 Kaheman ct al. 2020) introduces the concept of identifying implicit expressions of the form f(x. y) = in the context of differential equations (i.e. y = áे = da, for i E (1....,d)) d] where x E Rd). Further, Fureqa (Schmidt. & Lipson p009), well-cstablished bascline SR algorithm, focuses on discovering invariants rather then trying perform prediction directly. Inspired by the two aforementioned methods, and by the fact that the main objective of SR is to recognize correlations and define non-trivial interpretable models, GSR identifies rela- tions between the input and transformed output through scarching for expressions of the form g(v) f(x), which keeps the possibility open for predicting the output in a straightforward manner.
di
Computational complexity. Although genetic algorithms are inherently heuristic, understanding how our GSR algorithm operates and scales could still be valuable. Following Algorithm from Appendix A we can approximate the time complexity of GSR as:
ON Np (Ma ng= mm My ng + No + N + log Np))
(16)
where N, is the number of generations until the GP algorithın converges, and N₂ is the number of iterations until the ADMM algorithm converges. Recall that N, is the population size, Ms and Ma denote the number of gg x mB and gg x basis matrices applied to x and y, respectively, and is the number of paired training examples. More details about GSR's computational complexity can be found in Appendix A
11

---

Published in Transactions on Machine Learning Research (12/2022)
Compared to GSR, the special version s-GSR adopts vanilla SR (where ( is simply y) with the same GP algorithm and coefficient optimization process (through ADMM) as GSR. Hence, s-GSR's time complexity can be approximated as:
O(N.. Np (M& . "B mB + N8 + N + log N₂)
(17)
Although GSR's computational complexily contains an additional term of ON, N2 M "B*), )) the number of GP generations .N produced by GSR is often much less than that of s-G,, which explains the runtime advantage of GSR over s-GSR shown in Table g
GSR's expression ability. The term Generalized in GSR mainly stands for its ability to discover analytical mappings from the inpul space lo a transformed output space through expressions of the form g [x). This generalizes the classical SR task of identifying expressions of the form y= f() (i.e. the latter is simply a special case of GSR with ( = y). In addition, the term Generalized can denote the fact that we constrain the search space to generalized linear models, keeping in mind that the search space could be confined to other generalized spaces. Note that, by finding relations of the form g(y) = =xx); the expression ability of GSR could resemble that of classical SR tasks which search for relations y fc(x) where the composition function fa(:) is defined as fal.) = h o f(-) h(f(-) with h(-) 9 '-) if g(-) is invertable, or function class of similar expression ability 9 ()) if (() is not invertable. However, GSR takes advantage of the fact that the target y is a scalar, and hence, we can apply many basis functions to (through g(-)) without much increasing the complexity of the expression. In other words, we can avoid searching for functions equivalent to g (-) in a space that could grow exponentially with the dimension of the input feature vector by simply searching for their corresponding inverse transformations applied to the scalar target variable. This concept, which happens implicitly in our algorithm provides an odge for GSR over traditional SR mothods in terms of runtime, complexily, and smoothness of the search space. In short,, GSR discovers simplified expressions by reducing redundancies in the search space, which greatly saves the computational complexity of the search process. It is worth mentioning that, in principle, GSR's concept of fitting g(v) = f(x (instead of y f fx)) could be applied in conjunction with other classical SR methods; this may require some modifications to their marameter/coefficion
optimization process.
GSR: a simple yet promising algorithm. GSR combines features and benefits from the usually dis- parate fields of system identification and genetic programming. On the one hand, SINDy methods use some LASSO-like approaches (or sequential thresholded least squares) to conduct their sparse non-linear regression for finding solutions that take the form of a linear combination of basis functions. On the other hand, evolutionary algorithms are effective in finding basis functions that achieve optimal solution. In other words, GSR combines well established evolutionary methods with more classical system identification meth- ods. Although each of the algorithm components are relatively simple, the overall GSR algorithm achieves promising experimental performance, highlighting new insights, which can open up new research directions for future improvement.
6 Conclusion
We introduce GSR, a Generalized Symbolic Regression approach by modifying the formulation of the conven- tional SR optimization problem. In GSR, we identify mathematical relationships between the input features and some transformation of the target variable. That is, we infer the mapping from the feature space to a transformed target space, by scarching for expressions of the form g(v) = fxx instead of y = fxx. We con- fine our search space to a weighted sum of basis functions and use genetic programming with a matrix-based encoding scheme to extract their expressions. We perform several numerical experiments on well-known SR benchmark datasets and show that our GSR approach is competitive with strong SR benchmark mothods. We further highlight the strengths of GSR by introducing SymSet, a new SR benchmark set which is more challenging relative to the existing benchmarks. In principle, GSR's concept of fitting g(y) = fxx could be extended to existing SR methods and could boost their performance.
12

---

Published in Transactions on Machine Learning Research (12/2022)
References
Sungsoo Ahn, Junsu Kim, Hankook Lee, and Jinwoo Shin. Guiding deep molecular optimization with genetic exploration. Advances in neural information processing systems, 33:12008-12025 2020.
Ali Al-Roomi and Mohamed E El-Hawary. Universal functions originator. Applied Soft Computing, 94: 106417, 2020.
Thomas Bäck, David Fogel, and Zbigniow Michalewicz. Evolutionary computation 1: Basic algorithms and operalors. CRC press, 2018.
Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. Neural symbolic regression that scales. In International Conference on Machine Learning, 2021.
Stephen Boyd, Neal Parikh, and Eric Chu. Distributed oplimizalion and stolistical learning via the allernaling direction method of multipliers. Now Publishers Inc, 2011.
Steven Brunton, Joshua 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 (15):3932 3937, 2016.
Bogdan Burlacu, Gabriel Kronberger, and Michael Kommenda. Operon c++: An efficient genetic program- ming framework for symbolic regression. In Proceedingg of the 2020 Genelic and Evolutionary Compulation Conference Companion, GECCO '20, pp. 1562-1570, New York, NY, USA, 2020. Association for Com- puting Machinery. ISBN 9781450371278, doi: 011145/337772933396669
Chen Chen, Changtong Luo, and Zonglin Jiang. Elite bases regression: A real-time algorithm for symbolic regression. In 2017 13th International conference on natural compulation, fuzay systems and knowledge discovery (ICNC-FSKD), pp. 529 535. IEFF, 2017a.
Chen Chen, Changtong Luo, and Zonglin Jiang. Fast modeling methods for complex system with separable features. In 2017 10th International Symposium on Compulational Intelligence and Design (ISCTD), volume 1, pp. 201 204. IEEE, 2017b.
Chen Chen, Changtong Luo, and Zonglin Jiang, Block building programming for symbolic regression. Newrocomputing. 275:1973-1980, 2018a.
Chen Chen, Changtong Luo, and Zonglin Jiang. A multilevel block building algorithm for fast modeling generalized separable systems. Expert Systems with Applications, 109:25-34, 2018b.
Miles Cranmor. Pysr: Fast & parallelized symbolic regression in python/julia, September 2020.
Miles Cranmer, Alvaro Sonzaez- Peter Battaglia, Rui Xu, Kyle Crammer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. arXiv preprint arXiv:2006. 11287, 2020.
Stéphane d' Ascoli, regression for recurrence prediction. In International Conference on Machine Learning, pp. 4520 4520-4536. PMLR, 2022.
Kamienny, Guillaume Lample, and Francois Charton. Deep symbolic
Fabrício Olivetti de França. A greody scarch tree heuristic for symbolic regression. Information Sciences, 442:18 32, 2018.
Fabricio Olivetti de Franca and Guilherme Scidyo Imai Aldcia. inceraation-transffrmaatt cvolutionary algorithm for symbolic regression. Evolutionary Computation, pp. 1-25, 2020.
Anthony William Fairbank Edwards. Likelihood. CUP Archive, 1984.
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression. arXiv preprint Xiv: 08892, 2019.
13

---

Published in Transactions on Machine Learning Research (12/2022)
Kadierdan Kaheman, J Nathan Kutz, and Steven L Brunton. Sindy-pi: robust algorithm for parallel im- plicit sparse identification of nonlinear dynamics. Proceedings of the Royal Sociely A, 474462222222202202277 2020.
MeerreAAexannrr Kamienny, Stéphane d' Ascool,, Guillaume Lample, and François Charton. End-to-end sym- bolic regression with transformers. arXiv preprint arXiv:2204. 10532, 2022.
Samuel Kim, Peter Lu, Srijon Mukherjee, Michael Gilbert, Li Jing, Vladimir Čeperić, and Marin Soljačić. Integration of neural network-based symbolic regression in deep learning for scientific discovery. JEEE Transactions on Neural Nelwoorks and Learning Systems, 2020.
Michael Kommenda, Bogdan Burlacu, Gabriel Kronberger, and Michael Affenzeller. Parameter identifccation for symbolic regression using nonlinear least squares. Genetic Programming and Evolvable Machines, 21 (3):471-501, 2020.
John R Koza and John R Koza. Genetic programming: on the programming of computers by means of natural selection, volume 1. MIT press, 1992.
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabrício Olivetti de França, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason H Moore. Contemporary symbolic regression methods and their relative performance. arXiv preprint orXiv:2107.14551, 2021.
Mikel Landajuela Larma, Brenden K Petersen, Soo K Kim, Claudio P Santiago, Ruben Glatt, T Nathan Mundhenk, Jacob F Pettit, and Daniel M Faissol. Improving exploration in policy gradient search: Ap- plication to symbolic optimization. arXiv preprint rrXii:2107.09155 2021.
Peter M Lee. Bayesian statistics. Arnold Publication, 1997.
Thomas Leonard and John SJ Hsu. Bayesian methods: an analysis for statisticians and interdisciplinary researchers, volume 5. Cambridge University Press, 2001.
Qiang Lu, Jun Ren, and Zhiguang Wang. Using genetic programming with prior formula knowledge to solve symbolic regression problem. Computational intelligence and neuroscience, 2016, 2016.
Changtong Luo and Shao-Liang Zhang. Parse-matrix evolution for symbolic regression. Engineering Appli- cations of Artificial Intelligence, 25(6): 1193, 2012.
Changtong Luo, Chen Chen, and Zonglin Jiang. A divide and conquer method for symbolic regression. arXiv preprint Xiv: 08061, 2017.
Niall M Mangan, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Inferring biological networks by sparse identification of nonlinear dynamics. IEFE Transactions on Molecular, Biological and Multi-Scale Communications, 2(1):52- 2016.
Georg Martius and Christoph H. Lampert. Extrapolation and learning equations. arXiv:1610. 02995, 2016.
arviv preprint
Trent McConaghy. Ffx: Fast, scalable, deterministic symbolic regression technology. Genetic Programming Theory and Practice IX, pp. 235 260. Springer, 2011.
Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondicj Čertík, Sergey B Kirpichev, Matthew Rock- lin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. Sympy: symbolic computing in python. PeerJ Computer Science, 3:e103, 2017.
T Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P Santiago, Daniel M Faissol, and Brenden K Petersen. Symbolic regression via ncural-guided genetic programming population sceding. arXio preprint arXiv:2111.00053 2021.
John Ashworth Nelder and Robert WM Wodderburn. Generalized linear models. Journal of the Royal Statistical Society: Series A (General), 135(3):370-38) 384, 1972.
14

---

Published in Transactions on Machine Learning Research (12/2022)
Michael O'Neill and Conor Ryan. Grammatical evolution. IEEE Transactions Evolutionary Computation. 5(4):349-358) 2001.
Patryk Orzechowski, William La Cava, and Jason H Moore. Where are we now? large benchmark study of recent symbolic regression methods. In Proceedings of the Genetic and Evolutionary Computation Conference, pp. 1183-1190, 2018.
Yudi Pawitan. In all likelihood: statistical modelling and inference using likelihood. Oxford University Press, 2001.
Brenden K Petersen, Mikel Landajuela Larma, Terrell N. Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Representations, 2021.
Markus Quade, Markus Abel, Kamran Shafi, Robert K Niven, and Bernd R Noack. Prediction of dynamical systems by symbolic regression, Physical Review E, 94(1):012214, 2016.
Shahab Razavi and Eric R Gamazon. deural--eework-ddiretteddd genetic programmer for discovery of governing equations. XrXiv preprint arXiv:2203.00 08808, 2022.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. International Conference on Machine Learning, pp. 4442 4450. PMLR, 2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. science, 324 (5923):81 85, 2009.
Robert Tibshirani.
Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), (8(1):267-288) 1996.
Tony Tohme. The Bayesian validation metric: a framework for probabilistic model calibrulion and validation. PhD thesis, Massachusetts Institute of Technology, 2020.
Tony Tohme, Kevin Vanslette, and Kamal Youcef Toumi. A generalized bayesian approach to model cali- bration. Reliability Engineering & System Safety, 204:107141, 2020.
Tony Tohme, Kevin Vanslette, and Kamal Youmef Improving regression uncertainty estimation under statistical change. arXiv preprint arXiv:2109.082/ 2021.
Leonardo Trujillo, Luis Muñoz, Edgar Galván-López, and Sara Silva. neat genetic programming: Controlling bloat naturally. Information Sciences, 333:21 43, 2016.
Silviu-Marian Udreseu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 2616),,,,, 2020.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feymman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. arXiv preprint arXiv:2006. 10782, 2020.
Nguyen Quang Uy, Nguyen Xuan Hoai, Michael O'Neill, Robert McKay, and Edgar Galván-Toópez. Semantically-basee crossover in genetic programming: application to veal-valued symbolic regression. Ge- netic Programming and Evolvable Machines, 12(2):91-119, 2011.
Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. Symboliegpt: A generative transformer model for symbolic regression. arXiv preprint Xiv::22006. 14131, 2021.
Kevin Vanslette, Tony Tohme, and Kamal Youcef- Toumi. general model validation and testing tool. Reliability Engineering System Safely, 195:106684, 2020.
Marco Virgolin and Solon P Pissis. Symbolic regression is np-hard. arXiv preprint arXiv:2207.01011 2022.
15

---

Published in Transactions on Machine Learning Research (12/2022)
Marco Virgolin, Tanja Alderliesten, and Peter AN Bosman. Linear scaling with and within semantic backproopagation--bseed genetic programming for symbolic regression. In Proceedings of the genetic and evolutionary computation conference, pp. 1084 1092, 2019.
CJ Wild and GAF Seber. Nonlinear regression. New York: Wiley, 1989.
Hengzhe Zhang, Aimin Zhou, Hong Qian, and Hu Zhang. Ps-tree: A piecewise symbolic regression tree. Swarm and Evolutionary Computation, 71:101061, 2022a.
Michael Zhang, Samuel Kim, Peter Y Lu, and Marin Soljačić. Deep learning and symbolic regression for discovering parametric equations. arXiv preprint arXi::2277000528 2022b.
16

---

Published in Transactions on Machine Learning Research (12/2022)
A Implementation, Hyperparameters, and Additional Experiment Details
Implementation. Our GSR method discovers expressions of the form s(v) f(x) where R, x E R and where the search space for f()) and g()) is constrained to a weighted sum of Ma and M2. basis functions (namely of-) and ())), respectively. We use a matrix- based encoding scheme to represent o(-) and eb(-) using basis matrices B and B of sizes ng= Be and Bg x 1, respectively (where mg+e = n, +2 and no is defined in Section 3.3.2p. Hence, in addition to the number of basis functions Mg and My. the parameters "Bo; Ye; and BB allect the complexity of the evolved expressions, and hence can be controlled to confine the search space (although also affects the complexity of the expression, it is given by the problem and cannot be controlled). Although more than one basis matrix can lead to the same basis function (sce Remark 3.3 the search space of basis functions (mainly <()) is still huge in general, and thus, enumeraling all the possible basis functions is not practical. Hence, we will rely on genetic programming (GP) for effective search process. A pseudocode of our GP-based GSR algorithm is outlined in Algorithm
Algorithm 2: GP Procedure for GSR
Input: Np. Rp, Ma. Ms. Lx; Ly
Output: I
function SOIVEGSR(Np, nps Mg. Mas fix: Cu)
Initialize population:
I(k) // Initialize the generation (or iteration) counter
(2),, Ix, (k))
k 0:
for i= to N2 do
[kk + GRNVRATERANNOONNNNVIIAAAA
My. fx. Cy):
/- Each individual [k(k) contains two randonly generated sets of Mj and My basis matrices respectively
end
Evaluate each individual [k( with respect to the fitness function;
/* For each individual Z,(k), form the matrix A,kk), solve for the optimal coefficients vector w;( < (K) ")), then compute its fitness
I(k) + SORTED(I(K));
// in ascending order of fitness // Increment the generation (or iteration) counter
while Stopping Criterion not Sutisfied do
< ++1;
UPDATECRIVERION(:
/+ Start with strict stopping criterion (e.g. very lov error threshold) and slowly relax it (e.g- gradually increase the error threshold)
Lg L3 <- CHOOSESUBLIBAANY(G...
/- Choose sublibraries C fx and b S by allovable operations used with and y
respectively
[11:mp)(A) (A) I[i:mp] (k- 1);
/a The Up fittest individuals of the previous generation are copied to the current new one for i = + 1 to N, do
u - CENNEATERNDOONNTEERRLL..
if u = then
(k) + REPRODUCE([ (k), (g);
/* Crossover based on the surviving individuals
+/
else if & then
(k) < MUTATE([ (k), (g. LJ);
/* Mutation based on the surviving individuals
t/
else
(k) + GENRAVVRRNDOOSIOONDUUAAAA Ma, 22, Cg);
/* Randonly generate a completely new individual
*/
end
end
Evaluate each individual I,(k) with respect to the fitness function; I(k) < (ZK))
// in ascending order of fitness
end
I' - (k)); end function
// return the fittest individual
The main inputs to our Ga-basod algorithm are Ny. np, Mot My. fxs and by Recall that N, is the population size and np the number of surviving individuals per generation. Cx and Cg are the libraries of allowable transformations that be used with and y, respectively. These libraries form the first.
17

---

Published in Transactions on Machine Learning Research (12/2022)
two rows of mapping tables, e.g. Table 13 and are defined by the benchmark problem. Note that the division operator is not part of our GSR architecture. That is, the main arithmetic operations used by GSR are -, For example, for the Nguyen benchmark dataset, the library of allowable operations is Co f,= x, cos, sin, exp, In} as shown in Table 19 In this case, we define Cx = {1, cos, sin, exp, In} (resulting in the mapping Table 2 and by = (1, o exp, In]. Regarding the stopping criterion, common terminating conditions for GP include: i) a solution reaches minimum criterion (e.g. error threshold), ii) the algorithm reaches fixed number of generations (or iterations), iii) the algorithm generates a fixed number of individuals (or candidates expressions), iv) the algorithm reaches plateau such that new generations no longer improve results, v) combinations of the above conditions. In our case, the algorithm terminates when the solution hits a minimum root- squuue error (RMSE) threshold. To accelerate termination, we slowly relax the error threshold by gradually increasing it. To avoid reaching a plateau and since we are dealing with a small population size as shown in Table 9 we enhance diversity (in the basis functions) by producing completely new individuals with probability 1/2 per generation (while performing crossover and mutation with probability 1/4 each per generation). To speed up the search process, we employ sublibraries C3 C Lx and Cg S C,y of allowable transformations, used when generating completely new individuals (or completely new basis functions in the case of mutation). For x, we mainly rely on three sublibraries which are the most common: a polynomial sublibrary Cpoly; a trigonometric sublibrary Luig. and the original library Lx itself. For the Nguyen benchmark example above, Lpoly = = ""} and Lirig = {1, cos, sin}. Note that power operators such as o², 3 would be included in these sublibraries if they were part of the original library defined by the benchmark problem. The function CHOOSESUBLIRRARY For example, assuming k is the generation (or iteration) counter, if k S 1, 500, each cycle consists of 70 iterations broken into three stages, the first stage consists of 15 iterations and assigns 2 < fxx the second stage consists of 25 iterations and assigns f* + polly, and the third and final stage consists of the remaining 35 iterations and assigns 4 + fiig.. This cycle repeats until & = 1,500, after which the cycle's size becomes 1, 500 iterations broken into three equal stages (i.e. 500 iterations per sublibrary). For the cycle consists of 20 iterations, in which we equally alternate between the polynomial sublibrary poly and the original library C, itself (i.e. 10 iterations for each sublibrary). Indeed, the use of sublibraries is only possible when the correspouding operations are included in the original library defined by the benchmark problem (e.g. Neat-6 and Neat-8 cannot use trigonometric sublibraries since {cos, sin} are not included in their corresponding original libraries, as shown in Table 19p. In addition, it is up to the user to specify the cycle's size and how to alternate between sublibraries, or even decide whether to use sublibraries in the first place.
works according to some cycle.
Hyperparameters. Throughout our experiments, we adopt the following hyperparameter values. For GP, we use a population size Np = 30, and we allow for np = 10 surviving individuals per generation. We perform crossover with probability P. = and allow for only 2 parents to be involved the process (i.e. new individuals are formed by combing basis functions from two randomly chosen parent individuals). We apply mutation with probability P., & and allow for 3 basis functions (randomly selected from an individual) to be mutated (i.e. to be discarded and replaced by completely new basis functions). We generate (completely new) random individual with probability P, = 2. For ADMM, we use a regularizer À = 0.4, a penalty p= 0.1. The algorithm terminates when the e-norr of the difference between the weight vectors from two consecutive iterations falls below a threshold of 8 = 1 Regarding initial conditions, we use
I 11" (where a~* = denoles a normalized vector), 20 = = [1 [1.... 1J",, uo = = [0....0'T. V++.+
wo 2
For GSR, we allow for a maximum of Mo = 15 basis functions s(-) for each expression of s(-) (this is the maximum number since some of the Mo basis functions will be multiplied by 0, ie.. at most we get Mo nonzero coefficients multiplying the basis fcuntions). To avoid overfitting and overly complex expressions, we allow for a maximum of M = basis function v-()) for each expression of g(-) (in this case the maximum and minimum are both 1 and g(-) will consist of a single basis function). It is worth noting that we use Ma for Sy1Set-11. Each basis v-(-) will consist of single transformation ng= 1. Each basis Φ(-) will be a product of N, transformations, where N, is a random integer between 1 and 3, iee nB E 21). For each of these N₂ transformations, the maximum total multiplicity of all the independent variables (or features) in an argument is a random integer between 2 and 5, i.e. ny E (2, 4, 5J. GSR terminates when a candidate expression achieves a RMSE lower than a threshold with a starting value of e = 10- (recall that this threshold is slowly relaxed during the process, e.g. by progressively multiplying it by a factor of v10 for every 1, 500 iterations). All hyperparameter values are summarized in Table g
18

---

Published in Transactions on Machine Learning Research (12/2022)
Table Hyperparameter values for GSR for all experiments, unless otherwise specified.
Hyperparameter
Symbol
Value
GP Parameters
Population size
Nº
30
Number survivors per generation
np
10
Crossover probability
Po
1/4
Number parents involved in crossover
2
Mutation probability
P
1/4
Number of bases mutate
3
Randomly generated individual probability
Pr
1/2
ADMM Parameters
Regularization parameter
X
0.4
Penalty parameter
0.1
Tolerance on the solution error
8
10-5
Initial guesses
Wi, 20s U4
2 1,0
GSR Parameters
Maximum number of basis functions 6(-) for each expression f()) Maximum number of basis functions (5.)) for each expression of g())
Ms M ny
15 1 (2,3,4.5) (1,2,3) 1 10-6
Maximum total multiplicity of all features in an argument
Number of tranformations multiplied together per basis o(-)
Bg
Number of tranformation multiplied together per basis v-))
B
Tolerance on the solution error (RMSE)
Computational complexity. Although genetic algorithms are inherently heuristic, understanding how our above, we can approximato
GSR algorithm operates and scales could still be valuable. Following Algorithm the time complexity of GSR
o ONo- (Mo BB mB + M BB 1 N, + O (fitness) ) O (Np log N₂)
+ N. (IN,
((N,-mp)) -
mpp) (Po (crossover) + Pm O (muuation) + Pr (Mo Bg mB + MM Bg
+ No + O (fitness) ) + O LogplogNN)))
(18)
where N, the number of generations until the GP algorithm hits the tolerance c, and N₂ the number of iterations until the ADMM algorithm hits the tolerance Note that performing crossover mutation operations takes (11) time (i.e. constant amount. of time), and computing the fitness (which calculates RMSE on N paired training examples) takes O(N) time. Also note that Po, Pm, and P, are probabilities which can be treated as constants. Hence, GSR's time complexity reduces to:
ON Np Np (Ma up m m My Bg + N₂ + + log log
(19)
Compared GSR, the special version s-GSR adopts a vanilla SR (where g is simply y) with the same GP algorithm and coefficient optimization process (through ADMM) as GSR. Thus, s-GS'ss time complexity can be approximated as:
ON.. No (Mo Bg* mg= + N₂ + N + log No))
(20)
Although GSR's computational complexity contains an additional term of O (N. Np MM. ng+), the number of GP generations N. produced by GSR is often much less than that of G-R. which explains the runtime advantage of GSR over s-GSR shown in Table 12
19

---

Published in Transactions on Machine Learning Research (12/2022)
Additional experiment details. For all benchmark problems, we run GSR for multiple independent trials using different random seeds (following the experimental procedure in Petersen et al. [002]) Mundhenk et al. (2021)). Table I0 shows the recovery rates of GSR against liteoature-epporeed values from several algorithms on the Nguyen and Livermore individual benchmark problems. We first note that, due to the wide domain of sampled input points imposed by Livermore-1 (i.e. [-10, 10)), we observed some instabilities in the solution due to the presence of the exponential function, which we decided to exclude from the library of allowable operations during the search process for this benchmark. In what follows, we provide explanations for the results shown in Table 10 Note that Livermore-5 is difficult to recover by NGGPPS+SLP/ /HER and the remaining methods as it contains subtractions. Subtraction is more difficult than addition since it is not cumulative. This is not an issue for GSR since both additions and subtractions are equally recovered through the sign of the optimal coefficients multiplying the basis functions. Livermore-10 10 and Livermore-17 are more challenging than Nguyen- 10 since they require adding the same basis function many more times (which is apparent through the poor recovery rates of the different methods). Fortunately, this is also not a problem for GSR since it can be easily solved by finding the right coefficient multiplying the basis function. Livermore-18 18 is more challenging than Livermore-2 and Nguyen-5 since it requires recovering the constant 5 without constant optimizer (which can be recovered as =t=+ztat2) For GSR, this can be recovered by naturally solving for the real-valued coefficient. The problem of the different methods on Livermore-22 lies in the constant 0.5, which requires finding z+r compared to GSR which simply solves for the optimal parameter multiplying z², We observe that GSR performs poorly Livermore-9 and Livermore- 21 compared to NGGPPS+SLI /HER. This can be due to the choice of hyperparameters (e.g. Ma, "ge, and ne) as well as the GP-based search process. These two benchmarks require finding the first 9 and 8 powers of x simultancously, respectively, which can be difficult to achiove by GSR, espocially that we only consider Ma = 15 basis functions o(-) per expression of f(+), as mentioned earlier. Note that polynomials were not an issue for GSR up the 6th order (i.e. Nguyen-4). We also tried experimenting with a 7th order polynomial (i.c. y x + 2 + x5 + x + x² + +2 + x)) and GSR achieved 100% recovery rate. We started observing decline in the recovery rate when we added the 8th power of z. In other words, Livermore-21 (the 8th order polynomial) seems to be the limit that GSR can reach with polynomials while Livermore-9 (the gth order polynomial) becomes very difficult to recover. It is worth noting that if the libraries for the Livermore-9 and Livermore-21 problems contained the square and cube operators f (as is the case for the Jin benchmarks described in Table 19 then GSR would easily recover these two problems. Finally, GSR not able to recover Livermore-7 (y = sinh(a)) and Livermoro-8 (y = cosh(x)) since both benchmarks require finding the basis function e D which cannot be expressed using our current encoding scheme unless it is available a transformation by itself. That is, the exponential operator e* is not enough to recover 1 e- using our current encoding scheme. Had the negative exponential operator e been part of the library of allowable operations defined by Livermore-7 and Livermore-8, GSR would casily recover these two benchnarks. As we can see for SymSet-1 (u = s<in((<)- - +), we added the operator e to the library of allowable operations (see Table 2 which made GSR's mission much simpler and it was able to recover the corresponding ground truth expression as shown Table 26 It is worth mentioning that, although ground truth expressions are not expressible, GSR was naturally able to recover the best approximations possible for Livermore-7 and Livermore-8, which turned out to be their Taylor expansions around O. GSR's typical output expressions were as follows:
Livermore-7:
0.51655 = + x) 0.48165x(( x x)
+
-0.048738(z + x + z) x + +x) 0.0022335 =x(()([xxxx)
y RJ x + 0.166x + 0.008652.5
x
2 x +
(the first three terms of the Taylor series of sinh(x) around 0)
3!
5!
(1)
2Livermrre-7 and Livermore-8 can be expressed as sinh(x: =-- and cosh -- respectively.
20

---

Published in Transactions on Machine Learning Research (12/2022)
Livermoro-8:
y0.8616y = -0.2154 0.04295..4 + =)æ= 0.035877(z x) x -.2154 - 0.2154
0.0012368(z x (2 x x x x) 0.042956x(z + x)) 0.25898x(z) - 0.2154
y 1+ + 0.5z² 0.0416x
0.00144.
x2
1+
(the first four terms of the Taylor series of cosh(r) around 0)
4!
6!
2 co)h(x))
We next perform a runtime comparison between GSR and NGGPPS on the Nguyen benchmark problem set. We run each benchmark problem and report the runtimes in Table II We find that GSR exhibits faster runtime than NGGPPS, averaging 2.5 minutes per run on the Nguyen benchmarks compared 3.2 minutes for NGGPPS. It is worth noting that although GSR is, on average, faster than NGGPPS, it still exhibits slower runtime on some problems (e.g. Nguyen-6 and Nguyen-9 in Table [1]. This is due to the randomness of the search process as well as the use of sublibraries as mentioned earlier in the Appendix. Indeed, the runtime depends on the stopping criterion or condition. For example, one can shorten the runtime further if the interest is just in an approximation rather than an exact recovery. Our GSR method recovers exact expressions in the order of few minutes.
We further highlight the strengths of GSR on the new SymSet benchmark problem set, and show the benefits of scarching for expressions of the form g(v) f(x) instead of y f(x). Typical expressions, with exact symbolic equivalennc, recovered by GSR are shown in Table 26 The key feature of GSR lies in its ability to recover expressions of the form g(v) = f(x). To better highlight the benefits offered by this feature, we disable it by constraining the search space in GSR to expressions of the form = f() (which is the most critical ablation). We refer to this special version of GSR as s-GSR. Note that most of the SymSet expressions cannot. be exactly recovered by s-GSR (i.e. they can only be approximated). We compare the performance of GSR against s-GSR on the SymSet benchmarks in terms of accuracy and runtime (see Table [2]. The results clearly show that GSR is faster than s-GR,, averaging around 2 minutes per the SymSet benchmarks compared to 2.27 minutes for s-GSR (i.e. 11% runtime improvement). In addition, GSR is more accurate than s-GSR by two orders of magnitude. This is duc to the fact that GSR exactly recovers the SymSet expressions across most of the runs, while s-GRR only recovers approximations for most of these expressions. Tt. is worth mentioning that on SymSet-1, SymSet-4, SymSet-5, SymSet-10, and SymSet-12, we observe mean RMSE values of the same order of magnitude between GSR and s-GSR, since these expressions can be exactly recovered by simply learning expressions of the form y = f(x). As GSR has to perform search to discover that g(y) is simply y for these expressions, it exhibits slower runtime than s-GSR in recovering these expressions (see Table 12
In addition, we compare GSR against several strong SR methods with similar (or better) expression ability. In particular, we experiment on SymSet with NGGPPS, PSTree, PySR, and gplearn (see Table 12 GSR is more accurate than all these mothods by three orders of magnitude, which further demonstrates the advantage of our proposed approach. As for the runtime, PSTree is the fastest method, averaging around seconds per run on the SymSet expressions, while maintaining solid accuracies. This comes as no surprise given its state-of-the-art performance SRBench datasets (La Cava et al. 2021). It is worth mentioning that on SymSet-16 16, all the methods (i.e. NGGPPS, PSTree, PySR, and gplearn) exhibited some instabilities in the solution over all independent runs. Hence, we excluded SymSet-16 for these methods in Table 12
21

---

Published in Transactions on Machine Learning Research (12/2022)
Table 10: Recovery rate comparison of GSR against literature-reportei values from several algorithms
the Nguyen and Livermore benchmark problem sets over 25 independent runs. The ground truth expressions for these benchmarks are shown in Tables 19 and 20
Recovery Rate (%)
Benchmark
GSR
NGGPPS + SLP ALP/HER NGGPPS
GEGL
DSR
Nguyen-1
100
100
100
100
100
Nguyen-2
100
100
100
100
100
Nguyen-3
100
100
100
100
100
Nguyen-4
100
100
100
100
100
Nguyen-5
100
100
100
92
72
Nguyen-6
100
100
100
100
100
Nguyen-7
100
100
48
Nguyen-8
100
100
100
100
96
Nguyen-9
100
100
100
100
100
Nguyen-10
100
100
100
92
100
Nguyen-11
100
100
100
100
Nguyen-12"
100
4
12
0
Nguyen Average
100
92.00
92.33
86.00
83.58
Livermore-1 1
100
100
100
100
Livermore-2
100
100
44
Livermore-3
100
100
100
100
66
Livermore-1
100
100
100
100
Livermore-5
100
40
0
Livermore-6
100
100
88
64
Livermore-7
4
0
Livermore-8
0
0
Livermore-9
88
24
12
Livermore-10 10
100
8
24
0
Livermore-11
100
100
100
92
17
Livermore-12 12
100
100
100
100
61
Livermore-13
100
100
100
84
Livermore-14 14
100
100
100
100
Livermore-15 15
100
100
100
96
Livermore- 16
100
100
92
12
Livermore-17
36
Livermore-13 18
100
48
56
0
0
Livermore-111 19
100
100
Livermore-22 20
100
100
100
100
98
Livermore-21
88
64
Livermore-22
100
92
84
68
Livermore Average
85.45
77.45
71.09
56.36
30.41
All Average
90.59
82.59
78.59
66.82
49.18
22

---

Published in Transactions on Machine Learning Research (12/2022)
Table 11: Runtimes of GSR vs. NGGPPS on the Nguyen benchmarks. these benchmarks are shown in Table
The ground truth expressions for
Runtime (sec)
Benchmark
GSR
NGGPPS
Nguyen-1
18.66
27.05
Nguyen-2
25.96
59.79
Nguyen-3
38.77
151.06
Nguyen-4
63.82
268.88
Nguyen-5
447.01
501.65
Nguyen-6
465.79
43.96
Nguyen-7
33.35
752.32
Nguyen-8
93.89
123.21
Nguyen-9
391.79
31.17
Nguyen-10
68.05
103.72
Nguyen-11
38.86
66.50
Average
153.27
193.57
23

---

Published in Transactions on Machine Learning Research (12/2022)
Table 12: Average performance in mean RMSE and runtime, along with their standard errors, for GSR against s- s-GRR and several strong SR methods the SymSet benchmark problem sets over 25 independent runs. The ground truth expressions for these benchmarks are shown in Table 20
Mean RMSE
Runtime (sec)
Benchmark
GSR
S+ GSR
GSR
8-GSR
SymSet-1
2.52 10 + 9.29
10
3.69 1
+ 1.62
10
49.94 + 2.61
43.81 +: 1.95
SymSet-2
4.28
10
= 3.35
10
2.34 10
+ 1.21 x 10
75.24 2.83
91.95 4.26
SymSet-3
6.58
10
1.79
10
1.22 10
8.17 10
66.34
2.51
88.14 3.49
SymSet-4
10 1 5.19
10-4
1,37 10
1.46 x 10
428.23
405.34 6.04
SymSet-5
3.63 10 + 9.92
10-6
4,98 x 10 + 7.37 10
71.86 = 4.12
64.19 3.82
SymSet-6
4.38 x 10 + 8.09
10
8.12 x 10 + 1.93 x 10
76.59 - 3.72
94.18 + 4.16
SymSet-7
2,39 x 10 + 1.71
10
6,83 x 10 - 6.56 x 10
93.61 = 2.54
109.94 = 4.98
SymSet-8
6,83 x 10 + 1.95
10-5
7,88 x 10 + 3.88 10
68.07 + 2.47
90.61 + 3.06
SymSet-9
3.57 10 - 3.41
10
6.32 10 - 1.39 x 10
57.36 = 2.68
83.18 = 2.09
SymSet-10
2.12 10 1 4.43
10
3.24 10
7.58 10
393.62 5.47
379.56 6.86
SymSet-11
1.43 10
9.97
10
9.39 10
6.78 x 10
124.38
187.49 8.35
SymSet-12
10 + 5.37
10
7.09 10
2.88 x 10
78.53
2.97
SymSet-13
2,31 X 10 + 1.83
10
9.13 10 1 3.28 10
96.14
114.85 4.61
SymSet-14
2.18 10 + 9.47
10
6.14 10 + 7.78 x 10-3
112.32 + 4.37
137.61 + 4.53
SymSet-15
1,81 x 10 + 4.75
x 10-5
6,71 10 - 3.57 x 10
46.97 - 2.33
85.07 + 4.14
SymSet-16
7.24 x 10
- 3.58
10
9.86 X 10
- 2.19
10
98.37 = 3.71
126.16 - 5.94
SymsSet-17
2.57 10
1
7.03
x 10
3.41 X 10
1
4.54 X 10
116.78 1 4.49
143.31 T 6.28
Average
2.66 10
1
1.59
10
+
5.27 x 10
120.84
136.19
4.56
Benchmark
NGGPPS
PSTree
NGGPPS
PSTree
SymSet-1
3,66 10 + 7.81
10-3
7.92 10 + 2.23 x 10
175.32 + 1.54
42.98 + 3.23
SymSet-2
4.75 x 10
+ 5.92
10
1.96
10
2.72
10
171.46
+ 2.03
24.68 = 1.05
SymSet-3
1.34 10
1 x
3073
4.01 10
167.11
1.69
21.68
1.89
SymSet-4
1.76 10
= 9.89
10
103x
1.61 10
171.85
2.01
1.09
SymsSt-5
4.21 10
1 3.58
3.19x10-1
9.94
10
177.98
1.57
13.09
0.37
SymSet-6
2.08 10"
+ 4.80
10
1.4x 10
+ 3.72
10
179.34
1.75
13.28 0.36
SymSet-7
2.43 10
+ 9.62
10
5.25 10
+ 7.62
10
133.89
6.95
11.78 + 0.24
SymSet-8
3,93 x 10 + 5.08
10
4.02 x 10
+ 6.42
10
169.46 + 1.81
11.54 + 0.41
SymSet-9
1.34 10 +
3.55
10
1.23 10 + 1.90 10
175.04 + 1.04
12.23 + 0.29
SymSet-10
3.32 x 10 + 3.09
10
9.91 10 - 9.59 x 10
178.27 + 1.75
11.06 = 0.17
SymSet-11
2.49 10 -
2.71
10
1.1x 10 - 2.24 x 10
166.91 + 1.66
21.26 = 0.53
SymSct-12
8.76 10
1
6.89
10
8.32 10
T 1.39
10
178.56
1.53
11.65 0.32
SymSet-13
2.19 10
1.77
10
8.69 10 + 1.71
10
126.39
7.88
9.92 0.23
SymSet-14
4.93 10
+
2.85
10 18
8.2x 10
1 8.50
x 10
99.85
3.81
9.95 0.21
SymSet-15
2.86 x 10 17
+
4.49
10 18
6.22 10
+
1.28
x 10
94.42
2.84
9.93 + 0.18
SymSet-17
9.16 x -0
+
6.52 x 10-3
5.44
1 5.51
x 10-3
171.25
1.55
10.46 1 0.19
Average
4.65 10
L
1.24 10
4.57 10
T 7.98 10
158.57
1 2.59
16.23 T 0.67
Benchmark
PySR
gplearn
PySR
gplearn
SymSct-1
1.10 x 10 1 3.01 x 10
4.45 10
1 2.91 x 10
61.21 + 19.72
140.44 1: 7.27
SymSet-2
5,75 x 10 + 2.28 x 10
3.42 10
1 5.12
10
119.55 + 39.39
145.84 + 1.67
SymSet-3
1.71 10 + 2.02 10
1.72x
+ 9.41
10-3
12.62
+ 0.69
169.38 + 1.19
SymSet-4
6.14
6.12 10
x
+ 5.77
10
79.39 - 1.98
248.12 + 9.71
SymSet-5
5.91 10
1.28 10
2.
1 1.76
10
69.47 1 0.49
190.83 1.41
SymSet-6
6.57 10
x.9
10"
1.31
10
67.45 = 0.67
215.11 6.69
SymSet-7
5.22 10
x
8.93
1.50
10
117.47
169.94 1.91
SymSet-8
10 1
3.85 x
4.35
1 3.16 -
135.13 1 30.82
167.64 1 0.98
SymSet-9
2.58 10 +
8.51 x 10
2.19 10
+ 4.17
10
162.11 + 45.01
159.22 + 1.86
SymSet-10
4.1x 10 +
1.72
-
163 10
+ 3.38 x 10
48.17 + 0.43
192.92 + 1.72
SymSet-11
2,81 x 10 +
3.40 x 10
188 x 10
+
2.91 x 10
49.88 - 0.63
184.38 + 1.34
SymSet-12
2.53 x 10 +
1.11
X 10
2.62 -0
1
3.11
x 10
136.14 1 41.77
187.32 1 2.90
SymSet-13
8 + 6.39
10
1.07 10
+
1.29
10
16.01 - 2.98
12.52 + 0.77
SymSet-14
x.50 1
4.51
10
9.82 10
1.59
10
57.54 T 9.81
142.12
5.79
SymSet-15
6.90 10 +
4.02
10
10
1.36
x 10
93.54 =
147.58
SymSet-17
3.94 -
-0
1.18 10
8.20
10
167.45
65.81
148.42
0.67
Average
4.99 x 10 + 1.76 1
7.22 10
+
1.64 10
87.07 + 21.60
163.86 + 2.94
24

---

Published in Transactions on Machine Learning Research (12/2022)
B More Examples on Our Matrix-Based Encoding Scheme
The encoding process happens according to a table mapping rules that is very straightforward to under- stand and use. For example, consider the mapping rules shown Table I3 below, where the dimension of the input feature vector. Note that Table 13 involves more transformations than Table 2
Table 13: An example table of mapping rules for basis function. Placcholder operands are denoted by e.g. corresponds to the square operator. The identity operator is denoted by
by.1
2
3
4
7
Transformation (T)
sin
exp
lu
b..2
Argument Type (arg)
Σ
II
2
3
d
Variable (v)
skip
24 22
Example 1.
For d = 2, n= and mg+ (i.e. ne 2), the basis function 6(x) 29f²112 can be
generated according to the encoding steps shown in Table 14
Table 14: Encoding steps corresponding the basis function ((xx 272
Step
T
t1
Update
2
21
Tí(x) = of
2
exp
11
21
T2(x) = exiyy
Final Update: ((x) = (xx) T2(xx)
Based on the mapping rules in Table 13 and the encoding steps in Table 14 the basis function Φ(x) 2 can be encoded by 4 matrix as follows:
1 1 2
B
(21)
3, nB* = 5 and mg = (i.e. 3), the basis function ((x) 29 sin(araaa) (V
Example 2. For d =
can be generated according to the encoding steps shown in Table 15
201 +22
Table 15: Encoding steps corresponding to the basis function Φ(x) a5 sin( 2o2
20 2.14
Step
T
arg
Update
Σ
21
F2
Ti (x) = (2x) 182)
.2
T2(x) xx
J22
33
Tz(x) =in((122a)
v
Σ
22
13
Tx(x) = V222 + 210
Ts(x) =
Final Update: o(x)
T (x) Tx(x) Tx(x) (x))) Ts(x)
Based on the mapping rules in Table 13 and the encoding steps in Table 15 the basis function
o(x) 23 siu( alllaaaa\V//+22 can be encoded by a 5 matrix as follows:
22 +22
2
1
1
27
0
2
6
2
2
(22)
9
1
2
3
0
25

---

Published in Transactions on Machine Learning Research (12/2022)
Example 3. For nne 2, the basis function (b)) y² V" can be generated according to the encoding
steps shown in Table 16
Table 16: Encoding steps corresponding to the basis function vb4) - v²² VT.
Step
T
Update
Ti(y) 11a
2
T2(y) v
Final Update: <i(y)
Ti (v) Ta(v)
Based on the mapping rules in Table 13 and the encoding steps in Table 17 the basis function (0y) = y² V can be encoded by a x 1 matrix as follows:
[ B
(23)
9
Example For agv 3, the basis function vby) In(y) can be generated according to the encoding
steps shown in Table 17
Table 17: Encoding steps corresponding to the basis function (b()) In(y).
Step
T'
Update
1 Ti(y)
In Th(y) In(y)
1
Ta(y)
Final Update: <by) = (v) T2(y) Ta(y)
Based on the mapping rules in Table 13 and the encoding steps in Table 17 the basis function vhy)) In(y) can be encoded by x 1 matrix as follows:
0 B =
8
(24)
Example 5. For Agv = 1, the basis function (by) = et can be generated according to the encoding steps shown in Table 18
Table 18: Encoding steps corresponding to the basis function (b(y) = 6.
Step
T
Update
exp
Ty())
eV
Final Update: v(y) Ti(v)
Based on the mapping rules Table 13 and the encoding steps Table 18 the basis function (1()) = can be encoded by x 1 matrix as follows:
B [7]
(25)
26

---

Published in Transactions on Machine Learning Research (12/2022)
C Symbolic Regression Benchmark Problem Sets
Table 19: Specifications of the Symbolic Regression (SR) benchmark problems. Input variables are de- noted by x for I-dimensional problems, and by (2), for 2-dimensional problems. (a, indicates random points uniformly sampled between and b for every imput variable; different random seeds are used for the training and test sets. b,),b,c indicates evenly spaced points between and b for ev- ery imput variable; the same points are used for the training and test sets (except Neat-6, which uses E(1, 120, 120) lest set, and the Jin lests, which use U,-3,3, 3, as lest set). To simplify the notation, libraries (of allowable arithmetic operators and mathematical functions) are defined relative to 'base' library Co = -, cos, sin, exp, In). Placeholder operands are denoted by e.g. corresponds to the square operator.
Benchmark Expression
Dataset
Library
N
= 3 + x² +
U(-1,1,200)
Co
N
y x + 1 + x² +
U(-1,1,20)
Lo
Nguyen-3
y 2² + + + +2
U(-1,1,20)
Lo
Nguyen-4
y _6 + x + x22 +
U(-1,1,20)
Lo
Nguyen-5
sin(=")
U(-1,1,20)
Co
N
sin(x) + sin(x +
U(-1,1,20)
Lu
Nguyen-7
In(x + 1)+ In(2 +1)
U(0,2,20)
La
N
V
U(0,4,20)
Lo
Nguyen-9
sin(x1) + sin
U(0,1,20)
Lo
Nguyen-10 10
y (sin(rr)
U(0,1,20)
Lo
Nguyen-11
y x²²
U,0,1,20)
Lo
Nguyen- 12
= 2
x 1 12
U(0,1,20)
Lo
Nguyen- 12"
y 21
x
2
F2
U00,10,20)
Lu
Jin-1
2.5x
1.3r + 0.5x² 1.7x2
U(-3,3,100)
Lo [m
const
Jin-2
= 8x + 8n 15
U(-3,3,100)
Lo [In U
const
Jin-3
= 0.2xr + 0.5x 1.200 - 0.5x1
U( U(-3,3,100)
Lu {In] U
const
Jin-4
Ч 1.56. + cos((2))
U(-3,3,100)
Lo {In]
const
Jin-5
y 6sin(x)) c0)
3(-3,3,100)
Co {In]
const
Jin-6
y 1.35x + 5.5 ((x1 - 1) 1))
U(-3,3,100)
Lu [m]
const
Neat-1
= 4 + + +
U(-1,1,20)
61
Neat-2
= x + x² + 3 1² +
U(-1,1,20)
61
Ncat-3
У sin(x.) cos(2)
U(-1,1,20)
Neat-4
Ч In(x 11 In( +
U(0,2,20)
Neat-5
2 sin( cos(cc
U(-1,1,100)
Co
Neat-6
D 2
5(1,50,50)
Vol 1,,,,
Neat-7
2 2.1 cos (9.xz1 sin cos
E(-50,50,
Lo
U (tan, tanh, Vo}
-1)
Neat-8
U(0.3,4,100)
{+, exp, e
1.2-(2-2.5))
Neat-9
y 1+., 1+2,
B(-5,5,21)
Co
27

---

Published in Transactions on Machine Learning Research (12/2022)
Table 20: Specifications of the Symbolic Regression (SR) benchmark problems. Input variables are de- noted by x for I-dimensional problems, by (21, for 2-dimensional problems, and by T2, I3) for dimensional problems. U(a, b) indicates random points uniformly sampled between a and for every imput variable; different random seeds are used for the training and test sets. To simplify the notation, li- braries (of allowable arithmetic operators and mathematical functions) are defined relative to 'base' libraries Co = (+, X, cos, sin, exp, Inj or LE = Lo U [const]. Placeholder operands are denoted by ,, e.g. corresponds to the square operator.
Benchmark
Expression
Dataset
Library
Livermore-1
3 x + (x?)
U(-10, 10, 1000) U(-1,1,20)
Co
Livermore-2
sin(x²) cos(s)= ((((x)-) - 2
Lo
Livermore-3
sin(x" (z²) (0)(
U(-1,1,20)
Co
Livermore-4
In(x 1 In(x 2 1) + In(x)
U(0,2,20)
Lo
Livermore-5
#1 2 + T2
U(0,1,20) Lo
Livermore-6
4x 3x
+ 2n.²
U(-1,1,20) Lo
Livermore-7
sinh(x)
U(-1,1,20) Lo
Livermore-8
cosh(x)
U(-1,1,20) Lo
Livermore-9
_² + x²
+ x x x + x² x3
U(-1,1,20) Lo
Livermore-10
(i)( co8(22)
U(0,1,20) Lo
Livermore-11
U(-1,1,50) Co
21+-2
Livermore-12 12
U(-1,1,50) Lo
Livermore-13 13
U(0,4,20) Lo
Livermore-14 14 x² + +x² sin(x) + sin(x²)
U(-1,1,20) Lo
Livermore-15 15
.r
U(0,4,20) Lo
Livermore- 16
x²
U(0,4,20) Lo
Livermore- 17 4 sin (1) cos(22)
U(0,1,20) Lo
Livermore-15 18
sin (x²) cos(
U(-1,1,20) Lo
Livermore-19 19
_ + +xr + x² +
U(-1,1,20) Lo
Livermore-222 20
-==
U(-1,1,20) Lo
Livermore-22 21
_ + x
+ 8
x 33
U(-1,1,20) Lo
Livermore-22
0.5x²
U(-1,1,,00 Lo
SymSet-1
sinh(x))
U(-1,1,20)
2j (I)
U
SymSet-2
(x 3x 2.8z 5)
U(-1,1,20)
SymSet-3
(x 1.2,2 11.5)
U(-1,1,20)
28 (02
SymSet-4
0.8 -cos(x) + 4.2c sin
U(-3,3,20)
28
SymSet-5
4.5x7 + +x1zI 1.702 3.1
U(-1,,1,0)
C8
SymSet-6
U(-1,1,20)
321-2ž
SymSet-7
In + 4x
U(0,2,20)
28
SymSet-8
/ 5x + 14 2r +7
U(-1,1,20)
2 1
SymSet-9
(2xr +x2) 1
U(0,2,20)
LE
10
1.5 cos In(2122) - 2.5
U(0,1,20)
CE
SymSet-11
/2 cos 21) )= + 3 +4
U(-1,1,20)
{02)
SymSet-12
0.4rf + 6.2zz 3. 3.5.III3 - 4.5
U(-1,1,20)
2E
SymSet-13
2za
U(0,1,20)
x1t.3
SymSet-14
#13223 =1t
2(0,2,20)
SymSct-15
(x) + +222)
U(0,1,20)
2E
SymSct-16
U00,1,20
SymSet-17
in (0.2c tao
tm 0.5 cos (23)
U (0, 1,20)
C
28

---

Published in Transactions on Machine Learning Research (12/2022)
D Typical Recovered Expressions
Table Typical expressions (with exact symbolic equivalence) recovered by GSR for the Nguyen benchmark set. Note that the coefficients in the GSR expressions form unit vector due to the normalization constraint imposed by the Lasso problem in Equation 6 It is easy to verify (by simplification) that the GSR expressions are symbolically equivalent to the ground truth expressions.
Benchmark
Expression
Truth
y x² + x² +
Nguyen-1
(2.1395 X - +
-0.558
0.26972 000--0.1300 0.1395 xxxx)
GSR
0.2697x
0.2697x 0.0651 + x) 0.0651 1
0.1395(x x x x) 0.2697x 0.558 x
Truth
x² + x²3
x² +
Nguyen-2
0.58554y
GSR
-0.19518(z 2 x -
0.09759(z Z)æ 0.19518(x
Truth
y = x. x² x
x²
Nguyen-3
y.5 = 12. 0.25 x
x)æ
GSR
- 0.125
0.125(z x z) 0.125x
+
Truth
y = x + 1
-0.48318 -0.09663xr - - -. 0.16106 x x)
Nguyen-4
GSR
-0.24159 Z|. -0--0006660-00000000-1-
-0.24159 x 2 -
(2 0.096636z
0.096636r 0.096636r 0.16106(x x
Nguyen-5
Truth
y sin(x²) cx)
GSR
0.63246 0.63246 sin - 0.31623
Nguyen-6
Truth
y = sin(x) + sin( +
GSR
y = 0.5 0oao sin(
x))
Truth
y 1)
0.23652
Nguyen-7
GSR
10.17082z + 1 0117822 0 0.1095x(z + x) x)
+0277702Exx))-002222222777 x)x 0.23652
+0226811-444422-0078844440--212 +
Nguyen-8
Truth
y = V
GSR
0.83654 In(y) 111_00001111111111
Nguyen-9
Truth
y sin(x1) sin
GSR
0..77735 = 0 0.577335mmmm xx)
Nguyen-10
Truth
y = (s) cos(22)
GSR
0.44721 y 0.89 sin co)(
Nguyen-11
Truth
y = xx
GSR
0.70711 In(y) 0.7071
Truth
y = xf x 1 -2
N1
GSR
-0.4 -0.2(z2 x2) 0.43; 0.4z, + 0.4(2z) + mi) 0.4(z) x1
-0.401(m1 X x x1)
Truth GSR
== x
27 1/² 32
12
0.6g
0.1 1 22 1 x2) 0.3(m) + 21)(21 x x xx) + 0.3(z) x 21 x x
(2 x x1 0.622
29

---

Published in Transactions on Machine Learning Research (12/2022)
Table 22: Typical expressions (with exact symbolic equivalence) recovered by GSR for the Jin benchmark set. Note that the coefficients in the GSR expressions form unit. vector due to the normalization constraint imposed by the Lasso problem in Equation 6 It is easy to verify (by simplification) that the GSR expressions are symbolically equivalent to the ground truth expressions. Although GSR does not exactly recover Jin-6, it recovers approximations with very low RMSE, as shown in Table
Benchmark
Expression
Jin-1
Truth
2.6x 1.3x 0.523
1.72g
GSR
0.3U664y +0.153322
0.398632e]
0.26064drz 0.260614xz 0.76662,af
Jin-2
Truth
y = Sx} + 80.
GSR
0.06309 y -0.518175 0.1-22222220001111555+0002001110
0.276655
Jin-3
Truth
y 0.2x7 0.5xg 1.2xg
0.5z1
GSR
00.7943 0.11914m2
0.39775ur
Jin-4
Truth
1.5e + 5eos(z>)
GSR
0.25198y 0.37797e*
0.62995cos((22) 0.62995 cos(22)
Jin-5
Truth
y = 65in(22)
GSR
0.13484y
0.80904 sim( sin(2))
(i404522si(22 #1) 0.40152 +1)
Jin-6
Truth
1.35x + 5.5 ((rr 1)((22
1))
GSR
Not exactly recovered
Table 23: Typical expressions (with exact symbolic equivalence) recovered by GSR for the Neat benchmark set. Note that the coefficients in the GSR expressions form a unit vector due to the normalization constraint imposed by the Lasso problem in Equation It is easy verify (by simplification) that the GSR expressions are symbolically equivalent to the ground truth expressions. Although GSR does not exactly recover Neal-6, Neat-7, Neat-8, and Neat-9, it recovers approximations with very low RMSE, shown in Table 6
Benchmark
Expression
Truth
===" xf + + +
Neat-1
0.64952y = -0.32476xrz (z x 022000000000000000000000000--000
GSR
-0.32476(r x x) 0.2129x 0.32476(z 0117666
-0.064952( + 0.2129r O.27
Truth
y 2 + +
Neat-2
GSR
+0.18274z 0.18274x - 0.182744
Neat-3
Truth
sin(2 cos (os(s))
GSR
0.57735 0.57735 co)
0.57735
Neat-4
Truth
+ 1)
GSR
0.5 -0.55= x 0.5
0.25x 0.25%
0.5z( 0.25(z
Neat-5
Truth
y ===2: 2=)n((nnn co8(22)
GSR
0.57735 0.57735 cos(22) 0.57735
(.57733cos((c))siin
Neat-6
Truth
- E
GSR
Not exactly recovered
Neat-7
Truth
y = 2- 2.1 (en
GSR
Not exactly recovered
Neat-8
Truth
y 12- 2.577
-1
GSR
Not exactly recovered
Neat-9
Truth
1+3
GSR
Not exactly recovered
30

---

Published in Transactions on Machine Learning Research (12/2022)
Table 24: Typical expressions (with exact symbolic equivalence) recovered by GSR for the Livermore bench- mark set. Note that the coefficients in the GSR expressions form unit vector due to the normalization con- straint imposed by the Lasso problem in Equation 6 It is easy to verify (by simplification) that the GSR ex- pressions are symbolically equivalent to the ground truth expressions. Although GSR does not exactly recover Livermore-7 7 and Livermore-8, it naturally recovers their Taylor approximations, as discussed Appendix A
Benchmark
Expression
Livermore-1
Truth
+
sin(2²)
GSR
0.65079 0.21693 0.65079 sin(# 2) + x)
Livermore-2
Truth
sin(=")
GSR
sin)
0.166
Livermore-3
Truth
sin(x cos(
GSR
0.57735 y 0.57735 -
Truth
In(x + In
Livermore-4
GSR
-0. 0.50998 =
- 1| 0.022177
Livermore-5
Truth
24 x . 222
GSR
221 x1)
Truth
Ar 3r 2x x
Livermore-6
-0.15763
-0.26677x 0
GSR
-0.03804(z + [UL-0000---000000000000000000000
-0.093216(# 0099700--00--038811--------00000
Livermore-7
Truth
sinh(())
GSR
Not exactly recovered
Livermore-8
Truth
cosh((z)
GSR
Not exactly recovered
Truth
y == x + xx + 5
GSR
-0.30767 0..07777 x 0000
Livermore-9
-0.30767( x -100
0030767
X ( -
Livermore-10
Truth
6sin(m)) (o2
GSR
0.22912
0.68826 cos (22) sin(x1) 0.68826 (x) (05
Livermore-11
Truth
=
GSR
-0.40825 In(y) = - 0.8165 In x 01 0.40825 (22 + x1)
31

---

Published in Transactions on Machine Learning Research (12/2022)
Table 25: Typical expressions (with exact symbolic equivalence) recovered by GSR for the Livermore bench- mark set (cont'd) d). Note that the coefficients in the GSR expressions form a unit vector due to the normal- ization constraint imposed by the Lasso problem in Equation 0 is easy verify (by simplificatio)) that the GSR expressions are symbolically equivalent the ground truth expressions.
Benchmark
Expression
Livermore112
Truth
g
GSR 0.16903 In(y) (.845151m(()) In
Livermore-D 13
Truth
GSR
0.97332 In
0.1222
Livermore-1e 14
Truth
+ sin
GSR
0.40825 0.405555000005 0. si40825 0.400m0())0000000000000
0.40255
0,40825 sin
2)) 0.40825z
Livermore-Be 15
Truth
y == x
GSR 0.99015 0:0000053-1000000000005333333333 0.00000055|
Livermore-16 16
Truth
g =
GSR 0.99504 In (
Livermore-11 17
Truth
y (sin(an cos(2z)
GSR
0. 0.17408 0.69632
Livermore-15 18
Truth
cos(r)
GSR
0.19245
0.96225
Truth
2 +
Livermore-15 19
GSR
0146202
0. 20.200313
0.15297(2
0.15297(1
Livermore-20 20
Truth
8 = ==²
GSR
0.89442 In(y)= 0.47742
Truth
2 +x +
Livermore-21
GSR
-0.22497
Z) x ||g
Livermore-22 Truth
0.52
GSR
0.8165 In(y) = 0. 0..002252z
x/ 0.40825 x z)
32

---

Published in Transactions on Machine Learning Research (12/2022)
Table 26: Typical expressions (with exact symbolic equivalence) recovered by GSR for the SymSet benchmark set. Note that the coefficients in the GSR expressions form a unit vector due to the normalization constraint imposed by the Lasso problem in Equation It is easy verify (by simplification) that the GSR expressions are symbolically equivalent to the ground truth expressions.
Benchmark
Expression
SymSet-1
Truth
si)h(())
GSR
0.70448
0.35224e"z
0.17612ce 0.563584ec 0.17612ze
Truth
(x 3x
2.8x +)
SymSet-2
-0.11335 = -0.23188r
0.11335
0z.0651 - 0.11335
0.23188r
GSR
-0.11335
0.11335 0.23188
0.11335(z
-0.11335 + 0.34005z(zz x))
Truth
= (xx 1.2n + 11.5)*
SymSet-3
GSR
0.1173
+00.65766
0.265762' 0.26677
0 0.44965
+0.02346(z + + + +x)=" 0.2345
0.44965
Symset-4
Truth
0.8 cos(2) + 4.2c sin(
GSR
0.22485 - 0.112425 (o) 0.94437 si x 0.122 001124 0.17988
SymsSt-55
Truth
= 4,52
+ 2 1.70z
GSR
- 0. 0.16964y (2 x
x -
0.76338
0.525884 07277888+0028888888 0.288388r2
SymsSt-6
Truth
5
321-2
GSR
0.84515
-0.10903 x2 22 03300 0.50709kz
SymSet-7
Truth
= In(ar + 4x)
GSR
0.482715 (. + 0.21883.2 0.482715(z) 6r
21) 0.218832
SymSet-8
Truth
/5r 14 2.0 +7
GSR
0.060634 0.30317(zf x x 0.8188702ff211 x xx22 + 0014444188000000333333333333
0.060334(222
SymSet-9
Truth
= (2m1+ xx2)
+
GSR
0.83205 In(y) 0.5547 In(a +
SymSet-10
Truth
1.5cos((1) lu( 2.5
GSR
0.2244yy-00.511 08666 In(21
cos(2))
SymSet-11
Truth
/ /2 [1 30
+
GSR
0.228568 0.028571 0..457366 0.85713e**
(0571442css22)
Truth
0.4z +
6.2m² 53 4.5
SymSet-12
GSR
0. 0.18888 y 0.18324 0.055152 (n)
x 0.57008(xs xx) 0.504928.² 0.18324
-U.18324
0.132 0.5099888
SymSet-13
Truth
y
252
GSR
0.57735 In(y) U.57735
0.5773 +x)
SymSet-14
Truth
GSR
0.57735
In(y) = 0.57735 02 x mn)
0.57735 In(2 + 2s + 21)
SymSet-15
Truth
(2) +22
GSR
0. 0.70711 = 0.707
In +
SymSet-16
Truth
= _bdd)_b((c))998ccoof
GSR
0.098035 Iu(y) 0.960743 cos(ca) 0 In m2) + 0.25489Lm
SymSet-17
Truth
Iu =1 0.5 (zx3)
GSR
0.88045
0.17609e F1 0.440225 cos(13 x x3)
33

---

Published in Transactions on Machine Learning Research (12/2022)
E Limitations
Although GSR achieves great results whether by recovering exact expressions or approximations with low errors, it still has several limitations:
A Absence of division operations. The primary limiting factor to our GSR method is that it still cannot handle divisions. This is due to the way we define our encoding scheme. In this current version, we only consider a weighted sum of basis functions where the basis functions are a product of transformations; no divisions are involved. We can overcome this issue by modifying the encoding scheme to include divisions within the basis functions (e.g. a negative integer in the first column of the basis functions implies division by the corresponding transformation, i.c. using Table 13 a first-column entry of -8 encodes the division It). However, this will significantly increase the lotal number of possible combinations in which we can form basis matrices. Due to the lack of divisions in its current version, GSR suffers on some benchmarks such as Neat-6, Neat-8, Neat-9, Livermore-7, Livermore-8. It is worth noting that GSR can recover some divisions with the help of the In e- operators (see Livermore-11, Livermore-12, 12, Livermore-20, Livermore-22. SymSet-2, SymSet-6, SymSet-9, SymSet-13, and SymSet-14). This is only possible when the original function consists of only one term (not a sum of terms).
Composition of tranformations. Another limiting factor to the current version of GSR is that it cannot recover expressions containing composite functions, such as y = (oo(()) + In(x). In this example, the basis function gcos(s) cannot be recovered by GSR due to our encoding scheme. Again, if In(x) was not there, that is, if the function contained the first term only, i.e. y (cos(z), then GSR can handle the situation by recovering in(y) cos(z). The benefits of using g(v) fxx) can be clearly observed the SymSet benchmark problems (espocially SymSct-16, and SymSct-17).
Choice of hyperparameters and search process. Throughout our experiments, we have observed that, for some benchmarks (such as Jin-6 and Neat-7), althought they are expressible by GSR, they were not fully recovered. GSR only recovered approximations for these benchmarks with very low crrors. This can be explained by two reasons: i) The choice of hyperparameters affects the search process, ii) Our matrix- based GP search process may not be very effective on these benchmarks, given the complexity of their corresponding basis functions, and thus they may require a huge number of iterations to be recovered. That is, if we keep our code running for very long time, we may be able to recover these benchmarks. This can be verified by expanding Jin-6 and Neat-7 as follows:
Jin-6:
= 1.35x + 5.5 ((I) - () - 1))
= 1.352 + sin (I1 I1 - T2 + 1)
1.35x + 5.5| sin (-21 - x2) cos (x + 1) + cos (-z) x2) sin (x1 + 1)
= 1.35x - 5.5 cos sin(a1 + 222 cos(x1 +in3sin(1) sin 1 + x22) x2)
(x))
xx)
+ 5.5 cos cos(æ1 + 12) sin()122) +5.5 sin cos(01 + x2) (2)22)
(x()
(x)x)
Neat-7:
= - 2.1 cos(9.82)) sin (1.3x2)
2 - 2.1 (9.(9.8) cos(x1) - sin (9.8) sin(x1)) (sin(1.3) cos + cos (1.3) sin (222)
- 2.1 cos(9.8) s1.(1.3) cos(z,) cos(12) +2.1 (9.8) (1.3) sin cos(x2)
(x(x)
(x(x)
2.1 c9.(9.8) co.(1.3) cos(æ1) sin (x2) +2.1 (9. co3(1.3) sin (a sin(12)
(x)
(x
As we can see, GSR has to find the four corresponding basis functions simultaneously in order to recover the expressions.
Indeed, there are plenty of expressions that still cannot be fully recovered by our GSR method. This is the case for all the other methods well.
34


Reviews:
[{"review_id": "0hQPLi3Kwv", "reviewer": "Reviewer_zEUp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The proposed approach is well described in the paper. Its formulation is quite simple and reflects other previous works using matrices and genetic algorithms.\n\nThe definition of the method is well reported in the paper by some definitions. \n\nThe final part of the paper regarding the experimental evaluation should be extended considering other sota algorithms for solving the same problem. See for instance the paper \"Contemporary Symbolic Regression Methods and their Relative Performance\" for other methods and datasets.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents an approach called Generalized Symbolic Regression (GSR) trying to generalize the classical Symbolic Regression task. In particular the difference with respect to classical approaches is that the mapping is not directly done to the target variable but to a its transformation. \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "A more accurate experimental evaluation considering other methods and datasets could  improve the quality and the novelty of the paper.\n\nThe results are promising. As regards the time complexity it should be reported (sometimes the execution time could depend on the typical implementation and on the adopted programming language).\n\nAnother approach that can be used in the experimental evaluation should be that of reporting the complexity of the symbolic representation of the formula when compared to other methods.\n\nFinally, instead of using the recovery rate, it should be interesting to consider the symbolic solution definition reported in [1] (def 4.1) designed to capture models that differ from the true model by a constant or scalar.\n\n[1] Contemporary Symbolic Regression Methods and their Relative Performance", "broader_impact_concerns": "No ethical implications", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "M4dRIKurtP", "reviewer": "Reviewer_ZWfR", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- This work is generally well-organized and easy to follow.\n- Symbolic regression is an important and challenging problem in machine learning. Some TMLR audiences could be interested in this work.\n- The proposed GSR approach is straightforward and has achieved promising experimental performance. \n\n**Weaknesses**\n\nThe TMLR guideline asks the reviewer to check the following:\n\n(a) Are the claims made in the submission supported by accurate, convincing and clear evidence? \n\n(b) Would some individuals in TMLR's audience be interested in the findings of this paper?\n\nBased on the strengths of this paper, I believe the requirement (b) is well satisfied, but some claims made in this submission are not well supported.\n\n\n**1. The Claimed Generalized SR**\n\nThis work claims the proposed GSR approach generalizes the traditional SR approach, which searches for the form g(y) = f(x) instead of y = f(x). However, the proposed g(y) = f(x) formulation is very similar to the simple composition function formulation such as y = h(f(x)). If g() is invertible, we could have h = g^-1. Let f_new = hf, y = f_new(x) is the standard SR. What is the advantage of the proposed g(y) = f(x) formulation over y = h(f(x))? \n\nEven when g() is not invertible, it seems that it is still straightforward to find a function class h() such that y = h(f(x)) has a similar expression ability with g(y) = f(x). As correctly pointed out in this work, when g() is not invertible, the g(y) = f(x) formulation could indeed lead to an undesirable inference approach (e.g., a root-finding algorithm is needed). \n\nWithout a thorough analysis and discussion, it is hard to evaluate the actual advantage and contribution of the proposed GSR approach.\n\n**2. Linear Combination and Expression Ability**\n\nThis work restricts the f(x) and g(y) as a linear weighted sum combination of basis functions, but does not discuss its expression ability (with other methods). For example, given the same basic arithmetic operations and mathematical function, the formulation g(y) = f(x) should have better expression ability than y = f(x) (e.g., the s-GSR approach in the appendix). Therefore, it is expected that GSR should have better performance, especially when the ground truth function can be exactly represented by g(y) = f(x) but not y = f(x). \n\nHowever, the crucial question would be the comparison with a more fair SR counterpart with similar expression ability. It is always possible to improve the expression ability of y = f(x), such as (1) using the composition formulation y = h(f(x)), (2) enriching the mathematical function, and (3) relaxing the linear combination constraint. Will GSR still outperform those SR counterparts? Does the current better performance of GSR actually come from its new formulation that cannot be achieved by the SR approach with similar expression ability?\n\nThe expression ability for other SR methods are not analyzed or discussed in this work.\n\n**3. Experiment Setting and Comparison with Other Methods**\n\nThis work claims the current SR benchmarks (e.g., Nguyen, Jin, and Neat) are too easy for GSR and cannot really reflect its strengths. In this case, why not validate GSR's performance on the more realistic and strong SRBench benchmarks[1]?\n\nGSR is compared with four SR methods, namely, NGGPPS, DSR, BSR, and Neat-GP. However, according to [1], the other SR methods such as Operon [2], AIFeynman [3], and AFP_FE [4] could have much better performances for different problems. Is there any reason to choose NGGPPS/DSR/BSR/Neat-GP as baselines, and not compare with Operon/AIFeynman/AFP_FE or other methods with strong performance in [1]? \n\n**Other Comments**\n\n- The proposed GSR has 17 different hyperparameters. Will its performance be sensitive to those hyperparameter settings?\n\n- The function expression found by different SR methods will have different formula complexity, even when they are exact for the same ground truth expression. Will the GSR formulation lead to a higher formula complexity than other methods? \n\n**Reference**\n\n[1] Contemporary Symbolic Regression Methods and their Relative Performance. NeurIPS 2021.\n\n[2] Parameter identification for symbolic regression using nonlinear least squares. Genetic Programming and Evolvable Machines 2019.\n\n[3] AIFeynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. NeurIPS 2020.\n\n[4] Distilling free-form natural laws from experimental data. Science 2009.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes a generalized symbolic regression (GSR) approach to learn the mathematical expression from a set of independent variables to a target variable for a given dataset. The traditional symbolic regression (SR) approach usually has the form y = f(x) where f is a learnable analytical function. The proposed GSR is to learn the form g(y) = f(x) where g is also a learnable analytical function. \n\nFor the actual algorithm, this work formulates both g and f as a linear weighted sum of different learnable basis functions with a matrix-based encoding scheme. A genetic programming method and an ADMM-based optimization method are applied to 1) find the set of basis functions and 2) determine the weights for each basis function, respectively. Experimental results show GSR can achieve promising performance on some SR benchmarks and a newly proposed challenging SymSet dataset.\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address the major concerns listed in the weaknesses. It is crucial to:\n\n1) Analyze and discuss GSR's expression ability, and show why it could be better than the SR counterpart with similar expression ability such as y = h(f(x))\n\n2) Conduct more solid and convincing experiments to validate the performance of GSR with strong SR baselines on more realistic problems.\n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "RunzpCW2IN", "reviewer": "Reviewer_QcjW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The paper is well-motivated, and includes solid background and related work.\n- The experiments are solid and thorough: good comparison methods and benchmarks used.\n- The proposed method provides a first answer to an important open question: How to combine evolutionary methods with more classical system identification (sparse non-linear regression) methods.\n- Even though each of the algorithm components are relatively simple, the results are quite strong, highlighting many possibilities for future improvement.\n\nWeaknesses:\n\n- Although sparse system identification methods are cited, their technical connection is not mentioned in the methods section. These methods also use some LASSO-like methods for finding solutions that are linear combinations of basis functions. As is, it is not clear from the paper that such techniques are already well-developed.\n- Similarly, the reformulated version of SR is related to the key method in Eureqa of discovering invariants instead of trying to do prediction directly. Do the authors agree that this connection is related enough to be worth mentioning?\n- The new SymSet benchmark is useful to show importance of the SR reformulation, but it's not clear from the main paper how this benchmark was designed, e.g., what properties it is intended to have that other benchmarks do not. The reader would also like to see evaluation of the comparison methods on this new benchmark.\n- There are no plots of the GP convergence process, so it is difficult for the reader to get a feel for how the method behaves. \n- Although the proposed method is very practical compared to other SR methods today, it still avoids the key challenge in SR: how to find deep compositions of functions. The paper acknowledges this limitation, and says it will be addressed in future work, but does not suggest how it could be overcome, or why this new approach is a particularly good platform for handling this.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper develops a new approach for symbolic regression. The approach combines advantages from the usually disparate fields of sparse system identification (SINDy) and genetic programming (GP). This involves a new matrix representation for basis functions. GP is used to generate sets of such basis functions which are then plugged into a classical algorithm for finding sparse linear combinations of basis functions. The paper also introduces a reformulation of the symbolic regression (SR) objective that makes search smoother under their method, and makes their search space more expressive. The paper introduces a new benchmark that highlights the effect of this reformulation. In detailed experiments, the new method performs well compared to several key alternative SR methods on popular benchmarks.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Changes critical to securing my recommendation for acceptance:\n- Explicit mention of connections to related LASSO-like methods (e.g., SINDy methods) and finding invariants with Eureqa.\n\nChanges not critical to securing my recommendation for acceptance, but that would strengthen the work:\n- Move the \"Making Predictions\" section of the Discussion to the methods section, since the reader wonders this immediately. I'd put it at the end of Section 3.1\n- SymSet should be described at least at a high-level in the main paper, i.e., what are it's key properties that other benchmarks don't have? and what was the process for coming up with the problems in the benchmark?\n- It would be great to see some plot of how the evolutionary algorithm converges for at least some problems (at least in the appendix). This would give a much clearer sense of what the complete method is actually doing. Ideally, you could even show details of how the set of basis functions is refined over iterations.\n- The reader immediately is interested in seeing results for the alternate methods (NGGPPS, DSR, ...) on SymSet.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Three knowledgeable reviewers agreed that the proposed GSR is novel and the claims in the paper (post rebuttal) are supported by evidence. Only some minor points needs to be addressed in the revision before publication:
    
  - **Disentangling contributions** I agree with reviewer QcjW in staying that, in principle, one could apply GSR in conjunction with other classical SR methods. I recommend authors to (briefly) discuss this aspect in the paper. Authors can also mention that "Generalized" here is to be intended as in "Generalized linear models" as also pointed out by reviewers ZWfR and QcjW.  
  - **Add a discussion on computational complexity** I agree with reviewer zEUp that a discussion about computational complexity of GSR and competitors is missing. Especially one comparing GSR w.r.t. vanilla SR with the same genetic algorithm. I also agree with authors that big-O notation for worst-case complexity can be too pessimistic for genetic algorithms. However, it can be still valuable to understand how a method operates and scales.


Decision: Accept (poster)