Reinforcement Symbolic Regression Machine
Authors: Yilong Xu
Yang Liu
Hao Sun
Venue: ICLR 2024 poster (2024)
Keywords: symbolic regression
reinforcement learning
equation discovery
Primary area: neurosymbolic & hybrid AI systems (physics-informed, logic & formal reasoning, etc.)
TL;DR: Proposed a novel Reinforcement Symbolic Regression Machine (RSRM) that masters the capability of uncovering complex math equations from only scarce data.

Abstract: In nature, the behavior of many complex systems can be described by parsimonious math equations. Symbolic Regression (SR) is defined as the task of automatically distilling equations from limited data. Keen efforts have been placed on tackling this issue and demonstrated success in SR. However, there still exist bottlenecks that current methods struggle to break, when the expressions we need to explore tend toward infinity and especially when the underlying math formula is intricate. To this end, we propose a novel Reinforcement Symbolic Regression Machine (RSRM) that masters the capability of uncovering complex math equations from only scarce data. The RSRM model is composed of three key modules: (1) a Monte Carlo tree search (MCTS) agent, designed for exploration, that explores optimal math expression trees consisting of pre-defined math operators and variables, (2) a Double Q-learning block,  designed for exploitation, that helps reduce the feasible search space of MCTS via properly understanding the distribution of reward, and (3) a modulated sub-tree discovery block that heuristically learns and defines new math operators to improve representation ability of math expression trees. Binding of these modules yields the SOTA performance of RSRM in SR as demonstrated by multiple benchmark datasets. The RSRM shows clear superiority over several representative baseline models.

Published as a conference paper at ICLR 2024
REINFORCEMENT SYMBOLIC REGRESSION MACHINE
Yilong Xu', Yang Liu?, Hao Sun!.*
' "aolling School of Artificial Intelligence, Renmin University of China, Beijing, China; Schooo of Engineering Science University of Chinese Academy of Sciences, Beijing, China;
Emails: xuyilcag88e ruc. edu. en; 1 siuyang22eucas.com ac. on
haosun@ruc edu, on
ABSTRACT
In nature, the behavior of many complex systems can be described by parsimonious math equations. Symbolic Regression (SR) is defined as the task of automatically distilling equations from limited data. Keen efforts have been placed on tackling this issue and demonstrated success in SR. However, there still exist bottlenecks that current methods struggle to break, when the expressions we need to explore tend toward infinity and especially when the underlying math formula is intricate. To this end, we propose novel Reinforcement Symbolic Regression Machine (RSRM) that masters the capability of uncovering complex math equations from only scarce data. The RSRM model is composed of three key modules: (1) a Monte Carlo tree search (MCTS) agent, designed for exploration, that explores optimal math expression trees consisting of pre-defined math operators and variables, (2) a Double Q-learning block, designed for exploitation, that helps reduce the feasible search space of MCTS via properly understanding the distribution of reward, and (3) modulated sub-tree discovery block that heuristically learns and defines new math operators to improve representation ability of math expression trees. Binding of these modules yiclds the SOTA performance of RSRM in SR as demonstrated by multiple benchmark datasets. The RSRM shows clear superiority over several representative baseline models.
INTRODUCTION
The pursuit of mathematical expressions through data represents a crucial undertaking in contempo- rary scientific research. The availability of quantitative mathematical expressions depict natural relationships enhances human comprehension and yields more precise insights. Parsing solutions offer superior interpretability and generalization compared to numerical solutions generated by neural networks. Additionally, simple expressions exhibit computational efficiency advantages over the latter. As a result, these techniques have found applications across diverse fields, e.g., discovering fundamental physical laws (Udrescu & Tegmark [0200] [iu & Tegmark (202T) or governing equations (Schmidt & Lipson 2009; Chen et al. 2021 Sun et al. 2023), modeling material constitutive relations (Wang et al. 2019) and TCP congestion control Sharan et al. (2022), among many others.
The process of fitting expressions in early years involves polynomial interpolation to derive an equation, followed by the appearance of the SINDy method (Kaiser et al. (018). which utilizes sparse regression to identify appropriate mathematical expressions based on a predefined library of candidate terms. These methods c.g., (Sun al., 2021 Chen et al. 2021 Champion 2019), effectively reduce the search space from an infinitely large set of possibilities to a limited fixed set of expressions, thereby narrowing down the search process. However, the applicability of this approach is limited, since the compositional structure of many equations cannot be predefined in advance. Therefore, there is need for more comprehensive methods to search for expressions.
The Equation Learner (Martius & Lampert 2016 Sahoo et al. 2018) model was then introduced as a novel method in symbolic learning, which incorporates symbolic operators as activation functions. This modification enabled the neural network to generate more precise and interpretable functional relationships, allowing for the discovery of intricate math expressions. However, given the compact structure of EQL, optimizing the sparse network to distill parsimonious equations is a key challenge.
Another approach involves generating optimal expression trees (Hopcroft et al. (006)), where internal nodes correspond to operators and each leaf node represents constant or variable, By recursively
"Corresponding author

---

Published as a conference paper at ICLR 2024
computing the expressions of the sub-trees, these expression trees can be transformed into math expressions, Initially, genetic programming (GP) (Schmidt & Lipson 2009 Augusto & Barbosa 2000) Gustafson et al. 2005) was employed to address these problems. Although GP showed promise, its sensitivity to parameter settings leads to instability. Deep learning methods emerged then to tackle the problem. SymbolicGPT (Valipour al.] (022T) utilizes generative model like GPT to create expression trees, while AlFeynman (Udrescu & Tegmark] (2020) uses neural networks to analyze the relationships and dependencies between variables and search for relevant expressions. Despite its ad-hoc characteristic, the AlFeynman (Udrescu et al.) 2020) method was further improved, offering faster and more precise expression search capabilities. Additionally, reinforcement learning (RL) (Sun et al. 2023) has been employed, which utilizes the Monte Carlo tree search method to explore and discover expressions, along with a module-transplant module that generates new expressions based on existing ones. Deep (RI.) methods, e.g. DSR (Petersen et al. 2019). utilize recurrent neural networks to learn expression features and generate probabilities. A policy gradient search algorithm samples the probabilities to generate a batch of expressions, which are subsequently evaluated for performance. Combining DSR and GP leads to a new model called NGGP (Mundhenk et al.) 2021a), which achieves better performance. Then uDSR (Landajuela et al. 2022), comprehensive framework that combines DSR, AlFeynman, LSPT (Large-scale pre-training), GP, and LM (Linear models) emerged to enhance the efficiency and accuracy of symbolic regression. Pre-trained generative
models (Holt et al., 2022) and end-to-end transformer modules (Kamienny et al. 2022
Li et al.
2022 also achieved satisfactory expression search results.
Nevertheless, the existing methods still struggle with generating lengthy and complex equations, and are faced with issues related overfitting, c.g., poor generalizability. To overcome these challenges, we propose model named Reinforcement Symbolic Regression Machine (RSRM) that masters the capability of uncovering complex math equations from only scarce data, composed of an RL-search agent, a GP-based expression tuning element, and modulated sub-tree discovery (MSDB) block. The RL-search agent is designed based on the synergy between Monte Carlo tree search (MCTS) [Coulom|2000 and double Q-learning (Hasseli) 2010) for enhanced exploration and exploitation. The GP learner is employed to fine-tune the generated expression trees (c.g., see the demonstration in (Mundhenk al., 2021a)), while the (MSDB) block heuristically learns and defines new math operators to improve the representation ability of math expression trees.
We would like to emphasize that MSDB addresses crucial observation that models often struggle to generate complete expressions but excel in capturing certain components. For instance, NGGP (Mundhenk et al.
2a may discover an expression like -4-x² و3 + cos(y) + 1, while the ground truth is x== x² 0.5y + I. Notably, it successfully recovers the simplified expression x - 0.5y² with the same distribution. To this end, MSDB offers a new alternative to simplify expressions by subtracting specific components in the context of a sub-tree, as exemplified by the subtraction of 24 3 in the aforementioned case. Such an MSDB module takes the divide-and-conquee concept and could significantly improve the overall search performance of the RSRM model.
The aforementioned aspects form the main contributions of this paper: Our proposed RSRM model offers a novel solution to the search for mathematical expressions. By incorporating double Q-learning into MCTS, we effectively balance exploration and exploitation of SR tasks. The proposed (MSDB) block can handle equations with symmetry (reducing the complexity). and assist in dealing with long equations by identifying common patterns and defining new math operators on the fly. As a result, the RSRM model demonstrates clear superiority over several baseline models, which surpasses that of the baseline models in terms of accuracy and generalization ability.
BACKGROUND
Genetic Programming: Genetic programming (Stephens) 2016 Koza 1994: Schmidt & Lipson 2009) is employed to iteratively improve expression trees in order to approximate the optimal expression tree. The mutation step in GP enables random mutations in the expression tree, while genetic recombination allows for the exchange of sub-trees between expression trees, leading to the creation of new expression trees based on the knowledge acquired from previous generations. This "genetic evolution" " process progressively yields highly favorable outcomes after a few generations.
Double Q-learning (Hasselt 2010) is reinforcement learning algorithm designed to overcome the overestimation bias issue in traditional Q-learning. The key idea behind Double Q-learinng is to use two sets of Q-values independently estimate the value of each action in
Double Q-Learning:
2

---

Published as a conference paper at ICLR 2024
given state. By using two separate Q-functions, Double Q-learning can mitigate the overestimation bias of traditional Q-learning and provide more accurate value estimations, leading to better policy learning and performance in various reinforcement learning tasks.
Monte Carlo Tree Search: MCTS (Coulom 2006) is a decision- making search algorithm that constructs search tree representing possible game states and associated values. It employs stochastic simulations to explore the tree and determine the value of each node. This algorithm gained promi- nence via its adoption by the AlphaZero team (Silver al. 2017). MCTS consists of four steps in each iteration: (1) selection, (2) expansion, (3) simulation, and (4) backpropagation. During selection, the best child node is chosen based on certain criteria. If an expandable node lacks children, it is extended by adding available children. The simulation step involves simulating the current state before selecting the next node, often using the Upper Confidence Bound for Trees (UCT) algorithm to calculate the selection probabilities, defined as UCT('') = Q(0") + (N))) ((u") means the average reward of the child node, (v) and (vv) represents the number of visits to the current node and its child node, respectively, and typically represents the sxploration-exploitaaioo tradeoff parameter. The first part of the equation makes the nodes with high reward visit more often, and the second part ensures that the nodes with fewer visits have higher probability of being selected. Finally, in the backpropagation step, the reward function evaluates child nodes, and their values are used to update the values of parent nodes in the tree. The theoretical analysis (e.g., convergence, guarantees)) of the UCT-based MCTS algorithm can be found in Shah et al. (2022).
Here,
3
METHOD
best expressiions
④ [(r)
11x)
The RSRM model consists of a three- step symbolic learning process: RL- based expression search, GP tuning, and MSDB. With these steps, our model ef- fectively learns and represents the re- lationship present in the data, facilitat- ing accurate and interpretable modeling. The schematic representation of RSRM depicted in Figure [ The full settings of our model are in Appendix A
GP
MSDB
RL
R
2 R
Data reepproessiigg
ca lculate reward
convert s(x) odd or even function
add form
(xt
(l) - f(x) a
(y 07
find best g(x)
update
form
exp(x) emp(r) +1
form discovery
exp(2))
esp(s), exp(p)) logls)
The RL search consists of a double Q-learning empowered MCTS agent, Here. MCTS is employed for explo- ration (global search) that aids in gen- crating unexplored expressions, while double Q-learning enables exploitation that captures the local distribution of adopt
Figure 1: RSRM comprises several steps. First, the data is input into the data pre-processing module to determine the parity. Each epoch, the process initiates with the utilization of the MCTS module to generate expressions (Step D). Subsequently, these expressions undergo the addition of forms and reward calculation (Step 2). The expressions
are then refined through GP (Step 3)) leading to more ac- curate expressions (Step ④ Following this, the equation forms are integrated and rewards are calculated (Step 2)) This iterative cycle, involving MCTS and GP, continues for several rounds, ultimately resulting in the renewal of forms based on the most proficient expressions (Step )))
equations. Additionally, method that involves visiting each child node specific number of times before activating double Q Q-learning. This ap- proach aims to avoid excessive reliance on historical information, mitigating the risk of overfitting and promoting more robust learning process. To address the challenge of lengthy and hard equations, we introduce interpolation method (c.g.. data pre-processing to identify whether the equation exhibits symmetry prior to each search, followed by a modulated sub-tree discovery block (MSDB). If symmetry is present, we pre-process the equation accordingly to simplify the subsequent search process. This approach effectively reduces the difficulty associated with specific equations. The MSDB examines whether the few expressions that perform well adhere to a specific form. This divide-and-conguee algorithm enables a step-by-step search for equations, facilitating the generation of long expressions.
3.1 EXPRESSION TREE
The objective of SR can be transformed into the generation of an optimal expression tree (Hoperoft et al. 2006), which represents mathematical expression. The expression tree consists of internal
3

---

Published as a conference paper at ICLR 2024
Selection
Expansion
Simulation
Backpropagation
5
= 120
Select Max Reward
X 0.132 0.33
0.100 0.122
a' a'=UCT(n))
R 0.121
R 0.132
children
0.11 0.192
0.132 0.133
Random Pick
Pohildren
10%
10 simulations
R
0.129
R= [0.213, 1.13...
Double Q learning
,',d'RR
train
calculate reward after applying the form
Figure 2: Schematic of the proposed RL search. MCTS selects functions based on the maximum reward, expands them using the results of double Q-learning, simulates node selection through the UCT function, randomly fills the current tree, and provides rewards double Q-learning to train. Once the generation is complete, the rewards are back- poocppropagtted to the parent node.
nodes that correspond to operators (c.g., +, X, =, lo,,,, exp, sin, cos) and leaf nodes that correspond to constants (e.g.. 1, 2) or variables (e.g., x). By recursively computing the expressions of the sub- trees, the expression tree can be transformed into mathematical expressions. The process of generating an expression tree follows a recursive method where operators are added until no more can be added. This approach simplifies the task of creating expressions focuses on constructing the expression tree, which can be easily generated using recursive techniques.
In contrast to previous methods, employ a hierarchical traversal strategy for generating expression trees. This is motivated by the Monte Carlo tree search algorithm, where conducting more searches on vertices that are filled earlier deemed more beneficial. In the context of constructing expression trees, this implies that higher-level nodes in the tree carry greater significance. Consequently, We use hierarchical construction method build the expression tree layer by layer, similar to the hierarchical traversal of trees.
3.2 REINFOR LEARNING GUIDED SEARCH
The search step relies on the double Q-learning and MCTS algorithms, which are shown in Figure 2 The specific algorithm is shown in Algorithm
The reward function used in our approach is based on the root mean square error (RMSE) and is designed evaluate the fit of the generated equations to the measured data. It promotes concise and accurate expressions by assigning higher rewards to shorter and more precise functions. Inspired by the SPI. approach (Sun al. 2023), the reward function is computed by:
Reward function:
R E D:)
where is discount factor promoting concise trees, and / is the number of nodes in the expression tree. y: and yJ the true and the predicted values generated by the MSDB with the output Reinforcement Learning Search of the ith data point, respectively. Using this reward function, encourages the discovery of equations that minimize the RMSE and favors shorter and more concise expressions, leading higher reward values for functions that provide better fits to the data.
approach
4

---

Published as a conference paper at ICLR 2024
Algorithm 1 Expression generation by RSRM
dataset Siate expression form JF Parameters: discount rate 7). UCT const e, minimum selected times no Outputs: best expression
Input:
Initiate S as top of MTCS
Selection:
children of S with maxium R
Greedy selection
take action a
Simulation:
S" S
repeat
children of is empty then Expand S'
end if
if 3r E children of S" + (N(x) < no) then a' 6 else a' + randomly choose child of S' by UCT end if
D. Select child with visit times < 7ap b. Select through UCT function
S' take action S" S², Fill up randomly S" double Q-learning +- 5", a', R of S"
P train double Q-learning by simulated reward
until S' is full
Expansion
children of S double Q-Icarning + Pohstdren
estimate initial possibility of each child
Back-propagaae
Back-propagate R of 5' based on J.
Our method employs greedy selection, similar to Sun et al. (2023). Instead of selecting the token with the highest UCT score, we choose the token that currently yields the best reward (Eq. [ This ensures the selection of tokens leading to expressions resembling the current best one, potentially resulting in improved expressions, yet, increasing the possibility overfitting. Note that UCT is employed during the MCTS simulation while the greedy selection of the maximum reward is applied to choose the optimal expression tree.
Greedy selection:
Simulated reward: At each token generation, the entire expression tree is randomly completed based on the current tree. The reward is then computed using the reward function and fed back to double Q-laariig for training. This approach avoids excessive rounds of learning at the top node and filters out irrelevant nodes initially.
After expression tree is built, need to fill the parameter (i.e., equation coefficients) placeholders in it. We treat each placeholder as an unknown variable, which is optimized maximize the reward. The BFGS (Roger Fletcher & Sons 2013) algorithm, available in the scipy (Virtanen et al. 2020) module in Python, is used for optimization. In contrast to the approach in DSR (Petersen al. 2019), we find that Gaussian random numbers with a unit mean and variance provide more effective initial values for optimization (see further information in Appendix Section C.6 where we test the performance of the model with different initial values).
Parameter optimization:
3.3 MODULATED SUB-TREE DISCOVERY
We incorporate three specific sub-tree expression forms to enhance the exploration and analysis of cquations, where A represents a fixed form and f(x) a learnable part, explained as follows:
A+ f(x): This search form focuses on identifying expressions of the form like e - and ex + æ. By recognizing this pattern, we can effectively explore and analyze equations that follow the structure of e f())
A f(r): In this search form, we obtain good expressions such as 1.57c and 1.56c + aiming to detect equations of the form e" x (x))
Af(r); The search form fxx is designed to recognize equations like (e*)2.5 and (e*)s,
indicating the presence of expressions in the form ((a))
f(x)
Our approach involves the establishment of these forms based on the initial token of the expression tree, because the root of an expression tree serves as a focal point, indicating the primary operation or function in the expression. Thus, separate the sub-tree forms based it. Specifically, if the first
5

---

Published as a conference paper at ICLR 2024
token corresponds to addition (-) or subtraction (-), the method proceeds to learn the generation of the left and right sides of the respective operators. Similarly, for tokens such as multiplication (x), division (+), or exponentiation (*,, similar procedure is followed. In the case of unary expressions, such as trigonometric functions (sin and cos), the MCTS and GP models effortlessly derive the complete expression. Therefore, while our method involves degree of empirical design in identifying the sub-tree expression forms, it possesses a universal nature.
The complete form-discovery algorithm, which outlines the procedure for selecting and generating the search form among the three options, is provided in Algorithm 2 and Appendix Figure ST
Algorithm 2 Search for the form of the expression through the generated expressions
Input: best expression set Sheet
Parameters: selection ratio k. expression percentage ratio ky. maximum select number N Output: the form of the expression JF I + --enghhh of (Sbeat) Sort Sbeat by R decent
for i in 1, 2...! do
if < .V and R(Seaa [i]) 2 k Rmaz then
If number of & exceeds or R is low, break out Occurrences of function in Split-by-addiiioonSSS.. (b)
D D Split(Sbest())
end if
end for Gn f- D with maximum number of occurrences. if =C $ Z Go A then F = (())
The form is A Z means integer set.
else EC Z + Go x C then J A (()
D The form is A x (x)
else
F= f(x) for G in D do
b The form is A + F(x)
if Occurrences of 9 2 2 kp then J == F
> Add g to A
end if
end for end if
Splitting by Addition: In this step, convert the formula, which is represented as a token set, into a string using a library like sympy (Meurer et al. 2017). Then we expand the expression into a sum of simpler expressions. Next, split the expanded expression into multiple simple expressions using sum or difference notation. In this way, we convert [+, x, 2, y. 2, log, x| to x x y + 2 -log(z), and then transforms it xy, log(2).
Once the expression is refined into its desired form, the subsequent search becomes more manageable. For instance, in the case of aiming to derive exp(x:?) + xx + 3 + 0.5 log(z), can break down the search. Initially, we generate exp(x²) +.. to identify the form exp(x") + f(r), and then extend this exp(x=) x + x fxx) simplifying the process of obtaining exp (x + xt + x + 0.5 logg).))
Inspired by the approach proposed by Udrescu & Tegmark (2020). we introduce data pre-processing module determine the potential parity of the underlying equation. The cubic splines (Catmull & Rom 1974) are applied for equation fitting, generating a function. Subsequently, this function is used to compute the relationship between y(-z) and ((z), enabling the determination of whether y(x) is an odd, even, or neither function, where y(z) is the relation of x and y in given data.
When the error (RMSE) between y(-z) and y(r) remains below constant Esym= the function is considered even with respect to z. Negative values of the independent variables are transformed to their absolute values, while retaining the dependent variable values. Further exploration conducted using the form of ý (1(x) + g(-x))/2 to make it discover specific forms.
Similarly, if the error between y(-x)) and y(x) is within the limit constant Esym the function is classified odd relative æ. Negative values of the independent variables are converted to absolute values, and the dependent variable values are inverted. The search continues employing the form of ÿ (g(x) (x))// to make it odd in discovering specific forms.
Once the expression is refined to its parity form, the difficulty of searching for the expression is reduced. If we want get cosh(z), we only need to generate exp(æ) after parity determination. Such a partitioning strategy has been used in the past, e.g., Petersen et al. (2019) using sub-trees as new tokens, Sun et al. (2023) using transplanted sub-trees, Udrescu & Tegmark (2020) using problem

---

Published as a conference paper at ICLR 2024
Table 1: Recover rate (%) of several difficult equations in symbolic regression: trigonometric functions and sum of multiple power functions with parameter 1/2 in Nguyen; power functions and trigonometric functions in Nguyen", trigonometric functions and hyperbolic function and functions with weird power in Livermore; rational functions in R and R (where 0, y 0 add to dataset).
BenchMark
Equation
Ours
SPL
uDSR
NGGP
DSR
N
sin (2n (os(())))))))
100
95
55
72
Nguyen-12 12
2 - x 0.5x2 + +2
100
28
30
21
0
Nguyen-2
0.4841 + 3.39u] + 2.12 1.7821
100
94
100
98
90
Nguyen-9
sin(1.5m) + (0.50.)))))
100
96
0
65
Livermore-3
sin( (onx))ooo(xff) (x
55
15
0
0
Livermore-7
sinh (21)
100
18
0
24
3
Livermore-16
"1 sin (1 (un4))cos(<1))
100
40
60
26
10
Livermore-18
100
80
59
33
0
AlFeynman-9
+ 2
67
0
8
7
0
AlFeynman-10
=3 +2)
15
0
R-10
(I) (1) - +1
49
0
17
2
0
R-20
(x) +1/ +1)
89
0
0
0
0
R-3
(ar +
+
0
0
0
decomposition. Our main innovation is the use of partial sub-trees separated according to the plus sign as part of the new expression.
Also, we conduct parity determination performance test to compare the efficiency and effectiveness of form discovery by AlFeynman and our method. The experiment setting and results are given in Appendix G It shows that method (based cubic splines) outperforms AlFeynman (based on MLP) in terms higher accuracy and smaller data requirements.
RESULTS
Basic Benchmark Resu
100
Recovery raze (%)
We test the performance
Gum
of our method on multiple
SPL
40
UDSR
different datasets and com- pare it with the following baseline models in symbolic learning: SPL (Sun et al.
20
DSR
GP
Nguyen
Nguyen!
AlFeynman
2023). DSR (Petersen al.
SRBench Result
100
2019). NGGP MMundhenk
Recovery rate (9%)
al. 202 uDSR (Lan-
uDSR
dajuela et al. 2022). DGSR
DSSR
(Holt et al. 2022). gplearn
40
(Stephens] (016). and AFP.
20
AFP-FE
FE (Schmidt Lipson, (2010).
Feynman
Strogatz
Alla
Strogatze
The description of each base- line along with parameter set-
Figure 3: Recover rates of benchmark datasets. a, Basic benchmarks (detailed results shown in Appendix C). b, SRbench dataset (detailed results found in Appendix D and Appendix Figure (S3), where the symbol # denotes the presence of noise with mean of 10-3 added to the target values and represents missing data in the literature.
ting is found in Appendix B
4.1 BASIC BENCHMARKS
To evaluate the efficiency of our model, we first utilize four basic benchmark datasets (see Appendix C for details): Nguyen (Uy et al. 2011), Nguyen" (McDermott et al.) 2012). (Mundhenk et al. 2021b), Livermore (Mundhenk al. (021b), and AlFeynman (Udrescu & Tegmark (020). Note that parameter optimization (c.g., calibration of the equation coefficients) is prohibited in this experiment except the Nyugen dataset.
7

---

Published as a conference paper at ICLR 2024
We employ the recovery rate the evaluation metric, which measures the number of times the correct expression is recovered across multiple independent repetitions of a test. Note that this metric ensures that the model's 's output exactly matches the target expression.
We summarize the comparison of recovery rates for several difficult expressions on the five benchmark datasets listed in T The results demonstrate that model performs well on these complex expressions. Furthermore, we compared the mean recovery rates of all equations on each benchmark (see Figure 3 and Appendix C). Our method outperforms other approaches, achieving the highest recovery rates for all benchmark expressions.
We also conducted experiment on the trade-off between accuracy and the number of evaluations in the Nyugen Benchmark. Details of this experiment are given in Appendix Section C.7
4.2 SRBENCH DATASET
We further tested our model's 's ability to learn more complex equations with noisy training data using the SRbench dataset La Cava et al. (202T), where parameter optimization allowed. This dataset comprises 252 datasets sourced from Romano et al. (202T). We specifically concentrated on 131 of these datasets equipped with ground truth equations, which were drawn from two primary sources, namely, AlFeynman (Udrescu & Tegmark 2020) and Strogatz (Strogatz 2014) datasets.
Our evaluation encompassed a spectrum of baseline models from SRBench, as well as more contem- porary approaches like uDSR (Landajuela et al. 2022) and DGSR (Holt et al., (022). The recovery rate results of each method for different testing datasets (e.g., all datasets, only Feynman, only Stro- gatz) and noise effect can be found in Figure Here, the results for the Alfeynman dataset in Figure 3b are inconsistent with Figure since parameter optimization is prohibited in a but engraved b. Notably, our MSDB module demonstrated proficiency in handling intricate equations, such as the challenging example- - (32rt=3x² 333 + x. (a4))/(52555)) which can be effectively discovered through the form (C*+)// (2z)=)/(())))))))))))))))) x f(x). Additionally, numerous equations, like /2/ 10/3x1 1 + 1022. were successfully identified by aggregating multiple smaller equations, feat achievable through the form Cr + fx) Consequently, our model exhibited proficiency in uncovering wide range of equations. The comparative analysis in Figure 3p clearly illustrates the superior performance of our approach in comparison with other baseline models.
We further tested our model by discovering a surrogate formula to approximate the cumulative density function (CDF) of a normal distribution, e.g., (1, = f 1 where and a denote the mean and standard deviation. Since this equation lacks an explicit elementary expression, finding a parsimonious equation for approximation based on a small training set is intractable. Nevertheless, our model shows a great capability of approximating the CFD with generalizability (see Appendix Section for more details).
V2ra exp [-(t 100)//////// -x)²/(20?)] dt,
4.3 FREE-FALLING BALLS DATASET
We conducted an experi- mental evaluation on the free-falling balls dataset to assess the paramet- ric learning capability of our model. The dataset consisted of experimen- tal data of balls dropped from a bridge, as de- scribed in de Silva et al. (2020). The dataset com- prised 20-30 observations of a ball throw height within the first 2 seconds, aiming to learn the equa- tion governing the ball's drop and predict the height between 2 and seconds. Since an exact solution for this dataset is not available, employed the mean squared error (MSE) our evaluation metric.
Table 2: MSE of free-falling balls dataset. Details of the equations generated by different models are shown in Appendix
BenchMark
Ours
Ours'
SPL
M-A
M-B
M-C
baseball
0.053
0.068
0.300
2.798
94.589
3.507
blue basketball
0.008
0.027
0.457
0.513
69.209
2.227
bowling ball
0.014
0.034
0.003
0.33
87.02
3.167
golf hall
0.006
0.041
0.009
0,214
86.093
1.684
green basketball
0.094
0.045
0.088
0.1
85.435
1.604
tennis ball
0.284
0.068
0.091
0.246
72.278
0.161
volleyball
0.033
0.025
0.111
0.574
80.965
0.76
whiffle ball
0.038
0.660
1.58
1.619
65.426
0.21
whiffle ball
0.041
0,068
0.099
0.628
58.533
0,966
yellow whiffle ball
1.277
1.080
0.428
17.341
44.984
2.57
orange whiffle ball
0.031
0.368
0.745
0.379
36.765
3.257
Average
0.173
0.242
0.356
2.24
71.02
1.828
8

---

Published as a conference paper at ICLR 2024
We consider two sets of RSRM models, the standard one and the named RSRM (denoted by
Ours* in Table 2 that fixes the expression form c₁22 + C32² + I2r + C1 + f(z). We compared these models with the baseline method SPL (Sun et al., (023). since other models tend to have large generalization errors due to the limited data points (20-30 per training set) in the falling balls benchmark given the fact that the exact solution is unknown. Three physics models derived from mathematical principles were selected as baseline models for this experiment, and the unknown constant coefficient values were estimated using POWELL 1964). The equations of the baseline models are presented as follows. M-A: h(t)= = c2t² + C2t Cxt Ca. M-B: h(t) = ci exp(c2t) + Ost + Ca. and M-C: h(t) = c1 log(coshi (c2t)) + C3-
The results (see Table 2 show that in most cases, the RSRM model performs better than SPL. The RSRM model can successfully find the equation of motion for uniformly accelerated linear motion (cyx² + C2.: + C3 + f()) and search for additional terms to minimize the training error. This leads to improved results compared to SPI.. However, there are cases where RSRM makes mistakes, such as obtaining expressions in the form of e₁ cos(z)² + c2 + f(x) when searching for the yellow whiffle ball. This increases the generalization error and reduces the overall effectiveness compared to SPL. Overall, RSRM outperforms SPL in physics equation discovery, demonstrating its effectiveness in solving parametric learning tasks on the free-falling balls dataset.
Table Ablation recovery rate (%) of the Livermore dataset. The specific recovery rates are further shown in Appendix H
5 ABLATION STUDY
We conducted ablation studies on the Livermore dataset, namely. ablations of the double Q-learning (M-A), the
Equation
Ours
M-A
M-B
M-C
M-D
M-E
sin( (cos(an)
100
100
100
100
sin (cos(a])
20
55
0
MCTS algorithm (M-B), the MSDB (M-C), the pre-processing step (M-D), and GP (M-E), respectively. We list some of expressions affected by the performance of the model in Table
sinh((z)
100
100
100
100
10
100
100
83
100
88
100
67
100
100
100
67
100
100
100
100
100
12
100
33
Average:
97.95
94.36
93.64
80.45
89.45
84,95
When the double Q-learning module is removed, Model A with only MCTS experiences a decrease in knowledge from previous iterations. This results in reduced search efficiency but increased diversity. As result, we observe decrease in performance for equations like Σ z", while equations like (i)( - 1 show improved performance. On the other hand, when the MCTS module is removed, Model B with pure double Q-learning tends to overfit more quickly. Consequently, struggles produce the most challenging equations, such as sin sin(æ;) (=+))co((( cos((a))
1. Similarly, the absence of the expression form search module in Model C limits its ability to discover complex (in(=))oo((())) cos((?))
1. Lastly, Model D, without the preprocessing module, suffers a significant reduction in its ability to search for odd and even functions like s1nh( The removal of the genetic algorithm (M-E) resulted in decreased efficiency 1/3
expressions with simple forms, such x 1/3 and sin(
across all expression searches. While simpler expressions such as 2 still performed adequately, the performance for complex expressions notably deteriorated.
These observations highlight the importance of all the modules in RSRM. Each module contributes to the overall performance and enables the model to tackle different types of equations effectively.
CONCLUSION
We have proposed a novel model RSRM that integrates RL techniques, GP, and a modulated sub-tree discovery block improve the search process for mathematical expressions. Our model outperforms the -tatc-of-the-art basclines in the context of accurately recovering the exact equations for various datasets, and demonstrates superior generalization capabilities. However, one limitation of our current model is the lack of flexibility in setting the expression form as it currently encompasses only three fixed types, which restricts its adaptability to different problem domains. We anticipate future advancements in more flexible methods, e.g., potentially incorporating neural networks to generate slots for SR, utilizing other Q-learning techniques such as prioritized experience replay (Schaul et al. 2015) to enhance the exploitation, etc. Furthermore, we believe that our approach has the potential to be extended to other domains, such as reinforcement learning control tasks. By applying our method to diverse areas, we aim enhance the performance and applicability of SR techniques.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
The work supported by the National Natural Science Foundation of China (No. 92270118), which is greatly acknowledged. Code and models of Reinforcement Symbolic Regression Machine(RSRM) are available hhtpps///iihub.ccom/ wwwph//lloooo/mtllleeeeeaeeeeeee http:/ghhwh.hcoooom///////oo//// sci-comput/RSRM
REFERENCES
Douglas Adriano Augusto and Helio JC Barbosa. Symbolic regression via genetic programming. In Proceedinggs Vol. 1. Sixth Brazilian Symposium on Neural Networks, pp. 173-178. IEEE, 2000.
Edwin Catmull and Raphael Rom. A class of local interpolating splines. In Computer aided geometric design, pp. 317-326. Elsevier, 1974.
Kathleen Champion. From data to dynamics: discovering governing equations from data. PhD thesis, 2019.
Zhao Chen, Yang Liu, and Hao Sun. Physics-informed learning of governing equations from scarce data. Nature Communications, 12(1):6136, 2021.
Rémi Coulom. Efficient selectivity and backup operators in monte-carlo tree search. International Conference on Computers and Games, 2006.
Brian M de Silva, David M Higdon, Steven L Brumton, and J Nathan Kutz. Discovery of physics from data: Universal laws and discrepancies. Frontiers in artificial intelligence, 3:25, 2020.
Steven Gustafson, Edmund K Burke, and Natalio Krasnogor. On improving genetic programming for symbolic regression. 2005 IEEE Congress on Evolutionary Computation, volume pp. 912-919 IEEE, 2005.
Hado Hasselt. Double q-learning. Advances in Neural Information Processing Systems, 23, 2010.
Samuel Holt, Zhaozhi Qian, and Mihacla van der Schaar. Deep generative symbolic regression. In International Conference on Learning Representations, 2022.
J. F. Hoperoft, R. Motwani, and J. D. Ullman. Automata theory, languages, and computation. Pearson Education, 2006.
Eurika Kaiser, Nathan Kutz, and Steven L Brunton. Sparse identification of nonlinear dynamics for model predictive control in the low-data limit. Proceedings of the Royal Society A, 474(2219): 20180335, 2018.
Pierre-Alexande Kamienny, Stéphane d' Ascoli Guillaume Lample, and François Charton. End-to- end symbolic regression with transformers. arXiv preprint arXiv:2204. 10532, 2022.
John R Koza. Genetic programming as means for programming computers by natural selection. Statistics and computing, 4:87-112, 1994.
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabricio de Franca, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason Moore. Contemporary symbolic regression methods and their relative performance. In J. Vanschoren and S. Yeung (eds.), ), Proceedings of the Neurul Information
Processing Systems Track Datasets and Benchmurks, volume 1. Curran, 2021. URL https
oct.aeddde..s..eeere.......eerrr http://pppprffeeeeaaaaaaeeer files/paper/ 2021 22017//heeoccerrreeaacaee4666620
Mikel Landajuela, Chak Shing Lee, Jiachen Yang, Ruben Glatt, Claudio P Santiago, Ignacio Aravena, Terrell Mundhenk, Garrett Mulcahy, and Brenden K Petersen. A unified framework for deep symbolic regression. S. Koyejo, S. Mohamed, Agarwal, Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 33985-33998, 2022.
Wenqiang Li, Weijun Li, Linjun Sun, Min Wu, Lina Yu, Jingyi Liu, Yanjie Li, and Songsong Tian. Transformer- based model for symbolic regression via joint supervised learning. In The Eleventh International Conference Learning Representations, 2022.
10

---

Published as a conference paper at ICLR 2024
Ziming Liu and Max Tegmark. Machine learning conservation laws from trajectories. Physical Review Letters, 126(18)188004, 2021.
Georg Martius and Christoph H Lampert. Extrapolation and learning equations. arXiv preprint arXiv: 12....09.99 2016.
James McDermott, David R White, Sean Luke, Luca Manzoni, Mauro Castelli, Leonardo Vanneschi, Wojciech Jaskowski, Krzysztof Krawiec, Robin Harper, Kenneth De Jong, et al. Genetic program- ming needs better benchmarks. In Proceedings of the 14th annual conference on Genetic and evolutionary computation, pp. 791-798, 2012.
Aaron Meurer, Christopher Smith, Mateusz Paprocki, Ondfej Čertik, Sergey Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. Sympy: symbolic computing in python. PeerJ Computer Science, 3:e103, 2017.
Nathan Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio Santiago, Daniel M Faissol, and Brenden K Petersen. Symbolic regression via neural-guided genetic programming population seeding, arXiv preprini urXiv:2111.0005: 2021a.
Terrell Mundhenk, Mikel Landajuela, Ruben Glatt, Claudio P Santiago, Brenden K Petersen, et al. Symbolic regression via deep reinforcement learning enhanced genetic programming seeding. Advances in Neural Information Processing Systems, 34:24912-24923, 2021
Brenden K Petersen, Mikel Landajuela. T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXiv preprint arXiv: (9::::99222 1912.04871, 2019.
MJD POWELL. An efficient method for finding the minimum of function of several variables without calculating derivatives. The computer journal, 7(2):155-162, 1964.
John Wiley Roger Fletcher and Sons. Practical methods of optimization. 2013.
Joseph D Romano, Trang T Le, William La Cava, John T Gregg, Daniel J Goldberg, Praneel Chakraborty, Natasha L Ray, Daniel Himmelstein, Weixuan Fu, and Jason Moore. PMLB v1.0: open source dataset collection for benchmarking machine learning methods. arXiv preprint arXiv:20122000058v. 2021.
Subham Sahoo, Christoph Lampert, and Georg Martius. Learning equations for extrapolation and control. In International Conference on Machine Learning, pp. 4442 4450. PMLR, 2018.
Tom Schaul, John Quan, Ioannis Antonoglou, and David Silver. Prioritized experience replay. arXiv preprint arXiv: 1.55522, 2015.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324(5923):81--2) 2009.
Michael D Schmidt and Hod Lipson. Age-fitness pareto optimization. In Proceedings of the 12th annual conference on Genetic and evolutionary computation, pp. 543-544. 2010.
Devavrat Shah, Qiaomin Xie, and Zhi Xu. Nonasymptotic analysis of monte carlo tree search. Operations research, (6):3234-3260, 2022.
SP Sharan, Wenqing Zheng, Kuo-Feng Hsu, Jiarong Xing, Ang Chen, and Zhangyang Wang. Sym- bolic distillation for learned TCP congestion control. Advances in Neural Information Processing Systems, 35: 15:10684-10695, 2022.
David Silver, Julian Schrittwieser Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al. Mastering the game of without human knowledge. Nature, 550(7676):354 359, 2017.
Trevor Stephens. Genetic programming in python, with scikit-learn inspired api: gplearn. 2016.
SH Strogatz. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (westview, boulder, co), 2014.

---

Published as a conference paper at ICLR 2024
Fangzheng Sun, Yang Liu, and Hao Sun. discovery. Proceedings of the 30th International Joint Conference on Artificial Intelligence, 2021.
spline learning for nonlinear dynamics
Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discover- ing governing equations via Monte Carlo search. International Conference Learning Representations, 2023.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16)eaayy631, 2020.
Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. Advances in Neural Information Processing Systems, 33:4860-4871, 2020.
Nguyen Quang Uy, Nguyen Xuan Hoai, Michael O'Neill, Robert I McKay, and Edgar Galván- López. Semantically- crossover in genetic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Machines, 12:91-119, 2011.
Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. Symboliegpt: A generative trans- former model for symbolic regression. arXiv preprint arXiv:2106. 14131, 2021.
Pauli Virtanen, Ralf Gommers, Travis Oliphant, Matt Haberland, Tyler Reddy, David Coumaaeau,, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, al. Scipy 1.0: fundamental algorithms for scientific computing in python. Nature Methods, 17(3):261-222, 2020.
Yiqun Wang, Nicholas Wagner, and James M Rondinelli. Symbolic regression in materials science. MRS Communications, 9(3):793-800) 805, 2019.
12

---

Published as a conference paper at ICLR 2024
APPENDIX
A MODEL SETTING
In this section, we give more details about the settings of our models.
A.1 HYPERPARAMETERS
The full set of hyperparameters be in Table ST
Table S1: Hyperparameters of our model
Name
Abbreviation
Value
RL parameters
Minimum expression lengths
Imn
4
Maximum expression lengths
Ima
35
Maximum number of parameters
Cmax
10
Length discount rate
0.99
Training rounds
4
50
UCT constant
c
v2
Minimum selected times
ng
3
Learning rate of double Q-learning
Ir
10-3
Genetic Programming parameters
GP rounds
tgp
30
GP population
Pgp
500
GP number of best expressions
20
GP Mate rate
Pmnat
0.5
GP Mutate rate
Pmuutat
0.5
MSDB parameters
Error of Symmetry
Ezgm
10
Selection ratio
kg
0.1
Expression percentage ratio
kp
0.1
Maximum select number
N
A.2 EXPRESSION CONSTRAINT
In this study, incorporate a prior constraint inspired by the DSR (Petersen al. 2019) method to effectively reduce the search space for expressions. The following constraints are applied:
Length constraint: The length of expressions is restricted within pre-defined minimum and maximum values. If the current length falls below the minimum threshold, variables (., y. C, etc. ) and parameters will not be generated. Conversely, if the current length, combined with the number of nodes be generated, reaches the maximum length, only these nodes will be considered.
Unary operator constraint: The direct successor node of unary operator should not be the inverse of that same operator. This constraint ensures that the generated expressions adhere to the intended structure and prevent redundant combinations.
The successor node of a trigonometric function node should not be another trigonometric function. This constraint prevents the generation of expression structures that lead to unnecessary complexity or redundancy.
Trigonometric function constraint:
Maximum parameter limit: A specified maximum number parameters is imposed to control the complexity of the expressions and prevent overfitting.
13

---

Published as a conference paper at ICLR 2024
g
D
x
lxg(x")
x
100
Sort by R
+
log(x²)
xlog(x)
x
5
RL/GP
log(x²)
(2g(2x²)
x log(2x²)
xx
F = x xlog(x)
S best
Figure S1: Schematic of the proposed Form-discovery. We first obtain the output of RL/GP and select the equations in which the loss is relatively small. After that, we separate these equations by plus and minus signs and count the number of times these sub-expressions occur overall, and use the smaller equations with more occurrences to cobble together to create new expression forms.
By applying these expression constraints, we aim to enhance the search efficiency and guide the generation of meaningful expressions that align with the desired properties of the target problem.
A.3 GENETIC PROGRAMMING
Following the generation of expressions by each reinforcement learning algorithm, we engage in the optimization of predetermined set of expressions. For this purpose, we employ a genetic algorithm, utilizing the DEAP library in Python. This algorithm initializes half of the population using the outcomes of the reinforcement learning process, while the remaining half is generated randomly. Subsequently, we retain the most promising expressions and subject them to further analysis using subtree analyzer. This process serves update and refine the expression form, enhancing the overall efficacy of our approach.
A.4 FORM DISCOVERY
Presented below is an elucidative diagram (refer to Figure ST discovery outlined in Algorithm [2] This visual aid serves to elucidate the sequential progression within the algorithm and delineates the roles and interpretations of variables such as 5 Sayms G. D, and more.
pertaining to the process of form
B BASELINE MODELS
In this section, we give more details about the settings of baselines.
SPL (Sun et al. 2023): Rooted in the use of MCTS, SPL employs various greedy strategies ensure efficient exploration. utilizes full expressions sub-trees maximize the of past information. While excelling in shorter expressions, it falls short in handling longer ones.
DSR (Petersen et al. 2019 DSR takes a gradient-based based RL approach along with a recurrent neural network (RNN) that generates a probability distribution over expressions. While effective for longer expressions, it may exhibit limitations in generalization ability.
NGGP (Mundhenk et al. 202 Building upon DSR, NGGP enhances its capabilities. Expressions sampled via probability distribution undergo further optimization using GP. The refined expressions then train the RNN with risk-seeking policy gradient.
14

---

Published as a conference paper at ICLR 2024
uDSR (Landajuela al. 2022): The uDSR amalgamates DSR, AlFeynman, LSPT (Large- scale pre-training), GP, and LM (Linear models). excels in discovering formulas with constants (favoring polynomial type expressions), but albeit at the cost of increased compu- tation time.
DGSR (Holt al. 2022): This model leverages pre-trained deep generative models to capture the inherent patterns and structures within equations. This pre-training phase establishes a robust foundation for the subsequent optimization steps conducted through genetic programming.
The gplearn offers an efficient and rapid GP-based SR implemen- tation. While proficient in speed, it may exhibit instability and poor scalability.
splearn(Stephens 2016):
AFP-FE (Schmidt & Lipson 2010): technique that combines two important factors, age, and fitness, enhance the performance of evolutionary algorithms. FE means Co-evolved Fitness Predictors.
Age-Fitness Pareto Optimization is optimization
The full set of hyperparameters can be seen below.
SPL: In line with the original paper, we maintain the same parameter settings for SPL. The discount rate is set to 1) 0.9999, and the candidate operators include addition (+), subtraction (-), multiplication ( (x), division (+), cosine (cos(.)) )), sine (sin(-)), exponential (exp(-)), natural logarithm (log(-)), and square root (v-). Other parameter values are as follows: Maximum Module Transplantation: 20, Episodes Between Module Transplantation:
50000, Maximum Tree Size: 50, and Maximum Augmented Grammars: 5.
DSR/NGGP/UDSE In our study, we adopt the standard parameter configurations as pro- vided in the publicly available implementation of Deep Symbolic Optimization (DSO). This approach entails adjusting two primary hyperparameters. The entropy coefficient is set XH = 0.05 and the risk factor is set e = 0.005. Candidate operators are the same as those employed in the SPL.. Additionally, NGGP incorporates other hyperparameters related to hybrid methods based on genetic programming. The specific values are listed in Table S2
Genetic Programming (GP): We employ the gplearn library for GP-based methods. The hyperparameters for genetic programming are identical to those presented in Table [S2
DGSR/AFP-FE. For both of these models, we exclusively utilized the results obtained from the srbench dataset (La Cava al., (02T), and the parameters were meticulously tested in accordance with the specifications provided the official srbench dataset and its associated article.
Table S2: Genetic Programming Hyperparameters on baselines
Name
Value
Rounds
20
Population
1000
Mate rate
0.5
Mutate rate
0.5
C BASIC BENCHMARK RESULT
C.I OVERALL RESULT
In this section, give more details about the Min Depth and Min Complexity of difficult equations in each benchmark in Table S3
In this table, the performance of uDSR appears to be subpar, even falling short of its pre- decessor NGGP This outcome can be attributed to two primary reasons. Firstly, our appli- cation of symbolic learning lacks parameter optimization, except for the Nyugen-c dataset, in which uDSR notably excels due to its compatibility with parameterized scenarios. Sec- ondly, uDSR exhibits stronger tendency towards generating polynomial functions, whereas
15

---

Published as a conference paper at ICLR 2024
tests predominantly involve substantial number of trigonometric and exponential functions. For instance. expressions such as 2 2 + cos(zs)) result in complex equations like
a1 [ --0.043721$ 5.8423 0.01 3 0.0038xf25 0.492x442 + x4(222 + sin(xs))....
C.2 NYUGEN BENCHMARK RESULT
Nyugen Benchmark is standard benchmark for symbolic learning with two independent variables and equations randomly sampled over range. And Nyugen a parametric version the Nguyen benchmark, allowing the of parametric optimization test equations with parameters. In this section, we provide additional details about the results obtained from the Nyugen and Nyugen Benchmark experiment.
By referring to Table SS] readers can obtain more detailed information about the performance of each model on each expression, their comparative analysis, and any other relevant insights derived from the experiment.
C.3 LIVERMORE BENCHMARK RESULT
LiverMore Benchmark contains challenging equations rarely encountered symbolic learning. including high exponentials, trigonometric functions, and complex polynomials. In this section, provide additional details about the results obtained from the LiverMore Benchmark experiment.
By referring to Table S5 readers obtain detailed information about the performance of each model on each expression, their comparative analysis, and any other relevant insights derived from the experiment.
C.4 R BENCHMARK RESULT
R Benchmark consists of three built-in rational equations with numerous polynomials as divisors and devisees, increasing the learning difficulty. In this section, provide additional details about the results obtained from the R Rational Benchmark experiment.
By referring to Table S6 readers can obtain more detailed information about the performance of each model on each expression, their comparative analysis, and any other relevant insights derived from experiment.
Table S3: Minimum Depth of expression and Minimum tokens of expression Minimum Complexity of several difficult equations in each benchmark.
BenchMark
Equation
Min Depth
Min Complexity
5
sin x cos((:))
12
Nguyen-12
x x² 0.5x² + 82
24
Nguyen-2
0.48x + 3. + 2. 2.12x 1.78.r
25
Nguyen-9e
sin(1.5x1) + (0.50.))))
11
LiverMore-3
(in sin())cco((z1) sinh (x1)
15
LiverMore-7
15
2 2/5
LiverMore-16
7
17
18verMore-18
sin sin(x7)cos(m) 5
19
AlFeynman-9
+ 22 2 cos( +
17
AlFeynman-10
20
R-1
-
21
R-2
(x
+1/ +1)
31
R-3
35
16

---

Published as a conference paper at ICLR 2024
Table S4: Average Recovery Rate (%) of the Nyugen Benchmark over 100 parallel runs
Name
Equation
Ours
SPL
uDSR
NGGP
DSR
GP
Nguyen-1
at + af +
100
100
100
100
100
99
Nguyen-2
N
of mn of 2i
100
100
100
100
100
34
Nguyen-1
xf +
+.
5
sin(ar (oon(ss))) (21) -1
100
95
45
80
72
12
6
sin(ar)) sin(ar af)
Nguyen-7
lon(on 1) lox(sf +) Vel sin(an) t sin(a)) sin(21 sin(en)ood(ps) (ons(as)
100
100
100
35
17
Nguyen-8
100
Nguyen.9 9
100
100
91
100
100
17
Nguyen-10
100
100
86
Nguyen-11
"1 of 0.503 4.
100
100
100
13
Nguyen-12
Average
100,00:0.00
93. 13.50=111.7
87, 87,5
91.75=133.0
83. 131.8_188.4 18.5
44,424:2211
Nguyen-1
3,39,] 2. 2,12 187
100
100
100
100
Nguyen-2" 2
0.48xf 3.39xl + 2.12 + 1.78
100
100
100
100
0
Nguyen-5"
sin (xf()aco))) 0.75
100
Nguyen-7"
log 1.4) log(sf + 1.3)
100
100
2
Nguyen-8
/1.232
100
100
100
100
56
Nguyen-ye
sin(1.52)) + sin(0.5o))
100
Nguyen-10" 10"
sin(1 (o0o__5230
100
100
100
Average
100.00:0000
69.57:#35.2
60.71:+3.22
19.12
70.43=35.7
154
Table S5: Average Recovery Rate (%) of the LiverMore Benchmark over 100 parallel runs
Name
Equation
Ours
SPL
uDSR
NGGP
DSR
GP
Livermore-1 I
1/3 an sin(r:)
100
94
100
100
67
100
Livermore22
(io(zf))co((zz) -
100
29
58
61
26
1
Livermore-3
55
0
0
Livermor-4
log(2: + 1) log mx) log(x:)
100
61
100
100
72
100
Livermore-5
100
100
100
100
55
100
Livermore
Axf Bur 22 + smb((r))
100
100
100
100
100
Liveroore7
100
18
24
0
D
Livemmor-88
coblta))
100
Livermore99
S 2] Gein (oos(1ng)) (18o))/( 2s) af/fi
100
21
100
99
18
D
Livermore-10 10
100
100
100
70
23
Livermore-11
100
100
100
78
95
Livermore-12 12
100
100
100
100
Livemoor-133
100
100
100
59
0
Livermore- 14
+ + an + sin(srs) + sin(s?)
100
100
100
100
91
100
Livermore-1
100
28
2
Livermore-16 16
100
60
26
0
Livermore-17
Asinfrry (oon(ing) sin( sin(=f))con(a) (co)(())) -5
100
100
100
100
84
Livemore- 18
100
59
33
37
0
Livemmore119
+ an + +.0 exp(-xr)
100
100
100
100
Liveroore 20
100
100
100
100
100
100
Livermore-21
Y zi expl 0.5er)
100
52
100
100
13
12
Livermone-22
100
100
100
82
100
Average
97.95 4.0
51.00+16.9
14.6
80.68 14.0
10.41+15.7
50.77:-20.4
Table S6: Average Recovery Rate (%) of the R Rational Benchmark over 100 parallel runs
Name
Equation
Ours
SPL
uDSR
NGGP
DSR
GP
R--1
[m) 1>] mi
5
0
82
15
R--2
(r) Sz 1)/(2
80
0
40
0
R1-3
(uf + +++2)/(fff +
+ +1)
100
0
100
0
R--1
[i)
+ (2)
48
2
R--2
(2)
3z 1)(22 +1)
0
R133
+1
91
0
3
Average
68.83:28899
0.00=0.0 16.50=+26.2
26.67=31.1
0.00+00.0
0.00+00.0
C.5 AIFEYNMAN BENCHMARK RESULT
AlFeynman Benchmark contains lots of equations with physical meaning (as part of SRBench (La Cava et al. 2021)), such as expressions for gravity, kinetic energy, and light intensity superposi-
17

---

Published as a conference paper at ICLR 2024
tion. In this section, provide additional details about the results obtained from the AlFeynman Benchmark experiment.
By referring to Table S7 readers obtain more detailed information about the performance of cach model on each expression, their comparative analysis, and any other relevant insights derived from the experiment.
The selection 12 Al Feynman equations was based on a stratified representation of difficulty levels. The chosen expressions include easy ones, such 22 and 2 medium complexity expressions like #122243 (14( and 22 + 33g + #5 medium-hard expressions such 0.5a. (x + =3 + 3122za
and 4142T33 da ); and hard expressions like (4--,,14(66--)) ( and 1 + 1-2122/3 providing comprehensive coverage of the AlFeynman benchmark.
The comprehensive coverage of the AlFeynman benchmark reflected in Figure 33 of the results section, showcasing the outcomes for all AT Feynman equations considered in the problem set. Notably, our model achieves a state-of-the-ar 80% recovery rate, indicating its proficiency in capturing the underlying mathematical structures across a diverse range of expressions.
Table S7: Average Recovery Rate (%) of the AlFeynman Benchmark over 100 parallel runs
Name
Equation
Ours
SPL
uDSR
NGGP
DSR
GP
AlFeynman-1
100
100
100
100
100
100
AlFeynmun-2 2
100
99
97
100
97
AlFeynman-3
100
100
100
100
100
100
AlFeyaman 4
iggiga sin((rg)
100
98
90
100
100
AlFeynman-5
100
100
100
100
100
AlFeynman-6 6
an + ang coa(
100
100
100
100
100
AlFeeymnun-7 7
100
100
AlFeyanan 8
(ay irgat sin()a))
100
100
100
100
100
100
AlFeynman-9
#1 + ag + Aarnng (ns(
67
0
7
0
AlFeynman
(ag + af)
AlFeynman-11 II
0
0
AlFeynman-12
0
0
Average
73,50::24.1
66.42=27.8 8
63.50 63.50=26.0
67,25+27.4
66.42=21.8
60.58+22.7
C.6 CONST-OPTIMIZATION EXPERIMENT
We conducted an experiment considering different initialization approaches for constant optimization. Specifically, we explored seven methods:
Case 1. Initializing constants with vector of ones.
Case 2. Initializing constants with vector of random uniform values between 0 and 1.
Case 3. Initializing constants with vector of random uniform values between 0.5 and 1.5.
Case 4. Initializing constants with vector of random Gaussian values with a mean of 0 and standard deviation of 1.
Case 5. Initializing constants with vector of random Gaussian values with a mean of and standard deviation of 1.
Case 6. Initializing constants with vector of random Gaussian values with a mean of and standard deviation of 0.5.
Case 7. Initializing constants with a vector of random Gaussian values with a mean of 0 and standard deviation of //
We evaluated the recovery rate of each expression using different constant initialization methods across diverse benchmarks, employing ranges of input values such as [0,1], [-1,1], [0,10], 10], [0,50], [-00,00], and data sizes of 20 or 500. Notably, Table demonstraass that the average recovery rates across all initializing methods remarkably close.
Examining the distribution of expression recovery rates below 50%, 10%, and 0%, it becomes apparent that the method employing vector of ones exhibits the highest percentage in Table S9 This observation indicates that the vector of ones initialization method has the highest number of expressions unable to converge across 100 parallel runs.
18

---

Published as a conference paper at ICLR 2024
Table S8: Average Recovery Rate (%) of the Const Optimization Benchmark over 100 parallel on each size and each range
Name
Equation
Cisc ]
Case 2
Case 3
Case 4I
Case 5
Case
Case
Keijoes1
0 0.3 sin 28r1)
33.6
53.8
53.6
33.7
33.7
33.6
34.3
Keijor-4
(30.00 mazs)/((a 10.0)x2)
35.7
Knijee-14
8.0/ "1 +=)
100.0
90.8
100.9
82.0
95.8
99.6
99.8
Kenjoer-5
Nguyen-1
3.39 2.1221 7 1.78m1
97.2
97.5
96.8
95.9
96.9
Ngaeee-2"
0. 0.48x 3.30x] [ 2.12x 1.7851
98.0
92.3
97.9
92.9
95.0
97.5
96.2
Nguyen-se
sin(zf) con(x1) + 0.75
100.0
100.0
100.0
100.0
100.0
100.0
Ngaeeen?
log(2z 1.4) log(s; 1.3)
100.0
53.4
78.4
50.6
56.0
70,4
97.1
Nguyen-
v1.2321
100.0
100.0
100.0
99.5
Nguyen.ge
sin(1.521) 5x.5))
35.3
47.9
43.4
42.6
43.8
45.6
46.3
Nguyen-10
(in(1. cos(0.02g)
Nguyen-11
2.701
99.2
35.2
98.5
79.3
89.3
95.8
100.0
Jin" -1
2.52 1.zuf of 0.527 1.701
98.3
93.5
9403
97.0
94.7
91.7
97.4
Hio
8.0x 8.0s + 15.0e1
Fin 33
0.7 1.70r
100.0
54.1
92.2
52.9
76.1
83.3
87.9
lin
1.exp((1) (.0)))(21)
lio"
6.0 coa(111
100.0
Fin -6
1.35a 5.5 win( 172)
100.0
100.0
100.9
100.0
100.0
100.0
100.0
Average
81.06
76.19
80.00
73,45
77.17
79.74
81.08
Table S9: Persentage (%) of expressions under certain recovery rate about the Const Optimization Benchmark over each initializing method
Recovery Rate
Case
Case
Case 3
Case 4
Case 5
Case 6
Case
<50%
18.04
25.26
18.04
26.29
19.07
17.53
17.01
<10%
15.98
11.34
13.92
11.86
11.86
11.34
12.37
0%
7.22
4.12
6.19
3.61
4.64
3.61
3.09
C.7 TRADE-OFF EXPERIMENT
We used five different configurations for symbolic regression on the Nyugen dataset and obtained five different sets of results. We calculated the curves of the average number of tests about the expression recovery rate for all data/ Nyugen-4/ Nyugen-5/ Nyugen- Nyugen-12 12 for the five different data with different configurations shown in Figure S2
trade-off experiment
100%
Recover-rate
60%
ALL
40%
Nyugen-4
Nyugen-5
20%
Nyugen-11
Nyugen-12
50
75
100
125
150
175
200
Numbers Evaluations (1000x)
Figure S2: Trade-off between accuracy and number of evaluations in Nyugen Benchmark.
The evaluation count the RMSM exhibits a discernible pattern, roughly falling within two ranges: the first encompasses around 80% of the recovery rate, while the second spans from 80% 100%.
19

---

Published as a conference paper at ICLR 2024
Table S10: Trade-off experiment: Average Evaluation Number /Average Recovery Rate (%) of the Nyugen Benchmark over 100 parallel runs
Name
epochs
epochs
25 epochs
35 epochs
epochs
ALL
15654/53
27543/77
41525/87 87
55395/98
58240/100
Nyugen-4
28389/5
41638/43
54151/75
58483/100 100
76216/ 100
Nyugen-5
25084/0
57033/33
77075/80
118661/88
149868 100
Nyugen-11
20932/41
36487 /62
67237/88 /88
77757 100
84902 100
Nyugen-12 12
26345/0
65122/13
123571/41
175970/89
207889/95
Feynman
Strogatz
RSRM
AlFeynman
RSRM
AlFoynman
AFP
AFP FE
DSR
DSR
AFP
AFP
gpleam
gpleam
GP-GOMEA
GP-GOMEA
ITEA
EPLEX
EPLEX
Operon
Operon
SBP-GP
SBP-GP
BSR
Target Noise
BSR
Target Noise
FEAT
0.001
FEAT
0.001
FFX
FFX
0.01
MRGP
MRGP
20
60
50
100
50
100
Solution Rate (%)
Solution Rate (%)
Solution Rate (%)
Figure S3: Result of SRBench with 10 parallel runs for each dataset.
Both ranges be approximated with linear functions; however, the second range displays notably steeper slope, indicating that achieving higher recovery rates beyond 80% becomes considerably more challenging.
We categorized RSRM into using 5, 15, 25, 35, 45 epochs and tested the average number of tests and the average recovery rate for each different environment and different data respectively in Table S10
D SRBENCH RESULT
We used the full set of SRBench for testing, and the overall results are as follows Figure S3 including both the AlFeynman dataset and the Strogatz dataset.
E FREE-FALLING BALLS DATASET RESULT
In this section, provide additional details about the results obtained from the Falling-Balls Balls dataset experiment. To improve the performance of the SPL model, we incorporated the operators log(cor cosh(-)). This addition aimed to enhance the model's ability to capture the underlying patterns in the data.
The complete results of the experiment, including the functions found, can be found in Table S11
20

---

Published as a conference paper at ICLR 2024
Table S1 Functions generated in Palling-Balls Experiment
Name
Model
Equation
Ours
-4.43t + 0.36 + 1.51) + 47.35
baseball
Model A
0.09r
5.47t
2.47 46.52 + cos 2.5t)
SPL
0.6251 + 47.8 -4.44t
Ours
-1.66t* - 4.95/ cos -..95//900c0////////
46.46
blue basket ball
Model A
-0.10
4.49(2 37.54t 46.49 ((cos(1) 36.77)
SPI.
-0.25t* + - 46.17
Ours
-4.63/² sin
+ 46.13
bowling ball
Model A
0.18t 6.0r 6.0 2.12
0.62
SPL.
0.2856
3.82/ + 4.1 x 20000000+
12.45t + 46.1
Ours
-0.09r² - 4.44t + 5.26 10 log(t) + 49.51
golf ball
Model A
-2.18 + 11.75t + 1.96t 25.86 2.36 exp(t) 25.98 cos(t)
SPL
-4.9633/ + (og((cosh + 49.5087
Ours
46.34 4. 15t
green basket ball
Model A
-0.000 4.59t + 1.6t 45.26 +
020
+ 1) cos(t)
explcon(())
SPL
4.1465t + 45.9087 log(cosh(1)
Ours
47.78 cos c0. 0.02)
tennis ball
Model A
0.33t
1.9 0.661 17.74
SPI.
+ -4.05746
+ 17.8577
Ours
48.15 3.67(1 + 0.03)
volleyball
Model A
1.59r
11.1' + 0.931 58.53 c0.53
SPL
33.78 + 48.0744
Ours
(3.83 0.316) 47.07
whiffle balll
Model A
-0.00
2.17 1.69/ + 46.29 /L+ sin(3))
SPL
4.16/ + 47.0 exp( 0.15(")
Ours
2.181 + 0. (t)( 3.35 3.35 cus + 43.88
whiffle ball2
Model A
0.466
4.391 194 + 47.26 0.05 cos (exp((()
SPI.
65.86 (0.6xx((000076))) 18.61
Ours
(cos(1.75t) + 47.59) (0.
yellow whiffle ball
Model A
0.271
2.581 2.5t + 28 (1 0.11) exp( V + 2VP)
SPL
(148.99 14.58t" + 48.96 (og(cosh(c)))))))))o) (oggooo 3.065)
Ours
17..22 33.11/exp(()"-2 80.94
orange whiffle ball
Model A SPL
0.42t 3.811 -.44
1.666t 47..6
F GENERALIZATION EXPERIMENT RESULT
To compare the generalization ability of our model with other methods, we conducted an experiment on generalization performance. The dataset was generated using the cumulative distribution function (CDF) defined as
F(z, (L)) = 5 12
(1-a)) 2nd
dt
V 2TT
with varying means (μ) and variances (o). The dataset consisted of 201 points spanning the range from 1000 100. Each dataset is divided into three subsets: a training set, test set, and validation set. The training set comprises points ranging from 30 to 80, while the test set consists of points ranging from 10 25. The validation set covers broader range, spanning from 0 to 100.
Because this equation lacks an explicit elementary expression, it's impossible to obtain an analytical solution for the entire curve by reducing the training error to zero. Therefore, learning this function
21

---

Published as a conference paper at ICLR 2024
Curve-1 a = 1/30
Curve 2 o= 1/40
Curve-3 o= 1/50
1.0
1.0
1.0
0.6
extrapolatin
yd invoiivaaaatinnn
0.6
extrapelation
anni giviidaii
0.6
extrapolation
aini oiviidatio
Target
Target
Target
0.4
Durs
0.4
Durs
0.4
Durs
NGGP
NGGP
0.2
Cubic Spline
0.2
Cubic Spiine
0.2
Cubic Spiine
MLP
MLP
MLP
0.0
-100
-50
50
100
100
-50
0
50
100
100
-50
50
100
Figure S4: Result of generalization test experiment.
requires balancing between training error and generalization error, demanding a stronger ability to generalize. Moreover, only half of the curve data is provided, necessitating strong generalization skills from the model to extrapolate and accurately fit the missing portion of the curve. This calls for the proficiency of the model in capturing the distribution of the entire curve.
To evaluate the performance of our approach and compare it with baselines, we define the following settings for each method:
The training set is utilized for generating expressions and calculating the correspond- ing rewards. The test set is employed to evaluate the quality the generated expressions. Finally, the validation set is employed to select the most promising expressions from the outputs.
Ours:
(Gndhhunnheni al. 202 Both the training set and the test set are used for generating expressions and computing rewards. The HallOfFame, which contains the best expressions, is then leveraged to choose expressions using the validation set.
Linear regression: The training set and the test set employed for training the linear regression model.
Cubic splines:
The training set and the test set are used to train the cubic spline model.
Deep learning: the deep learning approach, we employ a Multilayer Perceptron (MLP) architecture with one input, one output, and a hidden layer ranging in size from 30 to 50. We set the learning rate 10-3 and train 100 epochs. We experiment with different configurations of the hidden layer and select the model that yields the best performance. The MLP trained using the training set, and the test set is used to evaluate the performance of each model. By varying the size of the hidden layer, we aim to find the optimal architecture that achieves the highest accuracy or lowest error on the given task.
The full results of the generalization experiment can be found in Table S12 and Figure S4 This table presents detailed overview of the performance of the model and the baselines. Additionally, Table presents the equations discovered by the model and the baselines. These tables demonstrate that our model outperforms the baseline methods in terms of generalization ability. The curves fitted by our model exhibit better accuracy and capture the underlying patterns in the data more effectively.
G PARITY DETERMINATION EXPERIMENT
In this section, we conduct an experiment to compare the efficiency and effectiveness of form discovery by
& Tegmark 2020) and our method.
According to AlFeynman (Udrescu Tegmark 2020), it also possible to use MI. as a learning curve determine the parity of function, compared the efficiency and effectiveness of MLP and cubic splines in determining parity.
The configuration of the MLP follows the structure employed in AlFeynman. It consists of a five-layer neural network with a balanced training set and validation set ratio of 5:5. The input layer accepts one variable and yields output of 128. Subsequent to the input layer, the hidden layers encompass input-output feature pairs 128-128, 128-64, and 64-64. The output layer produces a single variable. Optimization is executed using the Adam optimizer, and the activation function for each layer is set
22

---

Published as a conference paper at ICLR 2024
Table $12: Mean Squared Error (MSH) of each method and each part of the curve in the Generalization Experiment
Name
Ours
NGGP
Linear
Cubic Splines
MLP
total error curve 1
1.05 10
0.00215
0.0114
0.000381
0.0142
total curve 2
9.79 10
0.00163
0.0297
0.00162
0.261
total error on curve 3
2.61 10
0.327
0.0821
0.00563
0.0762
extrapolation
1.65 10
0.00429
0.0215
0.000758
0.0278
extrapolation error on curve 2
1.41 10
0.00324
0.0565
0.00323
0.518
extrapolation error on curve 3
2.67 10
0.65
0.158
0.0112
0.151
validation error on curve I
4.46 10-6
2.59 x 10 10
0.00129
2.76 x 10 10
0,000532
validation eτror curve 2
5.45 10-6
1.08 x 10 10-10
0.00265
2.87
10
0.00177
validation error curve 3
2.55 10-7
2.94 x 10-5
0.00518
3.75 10
0.00106
test error on curve test curve lest error on curve
6,95 10-6
9.09 X 10
0.0005455
< 1 10
0.000253
2.25 10-7
10
0.00127
< 10
0.00344
10-5
X
0.00272
10-12
0.00358
training error curve truining error on curve
3.88 10
3.57 10 10 12
0.000254
10 10
6.93 10-3
1.26 10-6
< x
0.000524
< 10
8.6
10
training error on curve 3
3.34 10-7
7.14 x 10
0.000905
< x 10-12
7.9 x 10-5
Table S13: Functions generated in Generalization Experiment. Our functions are easier to calculate and shorter than NGGP's.
Name
Model
Equation
Ours
0.503 + (117.088=) (/r 14702)
curve-1 1
NGGP
con(exp((0.4Be 20.50 + 2.78) 60.5))
Our
(6.08z 0.100/0000000
curve-2
NGGP wug.
Ous
curve-3
cos(log(1 16 4|
NGGP
0.6002)))
to the hyperbolic tangent (tanh). The training process encompasses 2000 rounds, initiated with learning rate of 0.01, which is subsequently reduced by a factor of 0. whenever the loss increases.
used seven functions, the top three are odd functions, from easy to hard, then the next three are functions, and the last one is non-odd non-even function.
We first tested the speed of both methods. The cubic spline method takes 0.216, 0.217, 0.219, 0.220, 0.226 milliscconds to Correspondingly. MLP takes 3.28, 3.64, 8.97, 64.1, 549.7 seconds to run.
10, 100, 1000, 10000, and 100000 data points, respectively.
We then tested two loss equations for the corresponding functions, shown in the Table S14 and the loss functions are as follows: Lodd = Din (v(x) + y((=))) Leven = Σ (x)x)= y( [x))²²²///
can be seen in Figure SS that 10 of MSE is used as the cutoff for whether it an odd/even function not, the amount of data required by MLP is about 100-1000 times more than that of the spline.
H ABLATION EXPERIMENT RESULT
In this section, we provide more detailed information about the results obtained from the ablation experiment LiverMore benchmark. The full results of the ablation experiment be found in Table (15) This table presents a comprehensive overview of the performance of the model under different ablation settings.
23

---

Published as a conference paper at ICLR 2024
Comparison of Parity Determination
10a
number of data need
spline MLP
10*
10a
10
odd-1
odd-2
odd-3
even-1
even-2 even-3
Figure S5: The parity determination performance test. We test the ability of parity determination of two models by three odd functions and even functions from easy to hard. Compete Setting of this experiment is shown in Table
Table S14: The average loss results for the Parity Determination Experiment across 10 parallel runs are presented in the table below. Each cell in the table contains two values: the upper value, situated above the horizontal line, signifies the loss value of the MLP, while the lower value indicates the loss of the cubic splines.
Name
Equation
Input Range
10
100
1000
10000
100000
[-,
26 10
10
1.15x10-2
5.88x10
1
odd-1
1.9x
10-1
2.78x10-17
3.47
1.95x 155x10-22
[-5,5]
3.0xx10
1.31x10-2
3.37x -
9.79x110
1.22
1.00100-11
3.18x10-1
0.00x10"
1.65x10 17
4.86x
[-,
2.15x1001
1,40x10"
1
odd-2
x + sinh(())) +
5,50010
[-5,5]
7.00100
4.5K100
117
2.96
2 +
[-1,
1,64x10-1
108x10-
5.76x10
2
odd-3
3.10X10-2
80
1.52x10
8.07x10-1
180x10-14
sin(x) cosh(a))
[-5,5]
7,60x10
7.07x10
1.95x10
155x10
3.45
816x10
1,28
275x10
182X10
[1,,]]
1
103x10
7.41x10
1,9x10
1676
even-1
2
10-
5
10
178x10
5,55
18
2.13x10
60
3
2 10
10-11
2.73x10-1
10
10
1.05X
[-1,]]
-
100
1035
2.65x10-4
7.32 10
even-2
x sinl((z)
10
12
2
[-5,5]
-00010
$.06110
288x10
3.15X
158810
1
1.00100
2.02
z log(z² + 1)+
[-1,1]
109
-
even-3
cos(x) x 10.
[-5,5]
9,99110
1000
9.21x10
6.65
7.76
1.20
161110
110X10
2.44100
1.09X TI
x log(2? 1)+ x7 + sinh(a)
[1,
9.90
108
5.18%
none
1.68
7.06x10P
1.23x
6.81x10
[-5,5]
4116110=
135x10
24

---

Published as a conference paper at ICLR 2024
Table S15: Average Recovery Rate (%) of the Ablation Experiment over 100 parallel runs
Name
Equation
Ours
ModelA
ModelB
ModelC
ModelD
Livermoe-1
1/3 + z1 + sin(x;)
100
100
100
100
100
Liverooe-2
sin(af (20)(2)) -2
100
100
100
100
Livermore-3
sin(ar (2)) (=r)) -1
20
55
Iog(21 +
log(sf 21) log(an)
100
100
Livermore-5
of
100
100
100
100
100
Livermore-6
1x
3x 2x +
100
100
100
100
100
Livermore77
sinh (n:
100
100
100
100
10
Livermore-8
cosh(ar)
100
Livermore99
83
100
Livermore-10
Gsin(ar (c0s(
100
100
100
100
Livermore-11 11
(#f=²) (un + ata)
100
91
100
100
100
Livermore-12
100
100
100
100
100
Livermore-13 13
1/3
100
100
100
100
Livermore-14
21 of 21 sin(a;) sin(2;)
100
100
100
100
100
1ivermore-15
100
100
Livermore-16
3/5 2
100
100
100
12
100
Livermore-17
4sin(ar) (20)
100
Livermore-18 18
100
89
100
Livermore-19 19
100
100
100
100
100
Livermore-20
expl
100
100
100
100
100
Livermore-21
100
100
100
100
100
Livermore-22
expl (l5z
100
100
100
100
100
Average
97.95:14.0
94.68
95.00:18.9
80.45:115.6
89.45:1: 11.9
25


Reviews:
[{"review_id": "cPTGJLnkWi", "reviewer": "Reviewer_jkfW", "summary": "This paper proposes Reinforcement Symbolic Regression Machine (RSRM), a symbolic regression method that combines reinforcement learning (RL), genetic programming (GP), and a novel modulated sub-tree discovery block. RSRM alternates between a reinforcement learning stage that searches for the optimal expression tree using double Q-learning and MCTS and a GP stage that refines the expression trees. At the end of each epoch, MSDS discovers new expression forms that reduce the search space for subsequent steps. RSRM demonstrates superior performance over baselines across a suite of symbolic regression datasets.", "questions": "- Can you elaborate on each step in figure 1? The way it's presented now makes it hard to understand.\n- What are some practical domains that SR is useful in? In particular, what are some scenarios where SR outperforms neural network regression?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The idea of combining RL and GP is a novel contribution. Using double Q-learning in junction with MCTS effectively reduces the search space. This is further enhanced by discovering new expression forms in MSDB.\n- The method achieves strong performance across a suite of datasets, outperforming baselines in terms of expression recovery rate.\n- The authors conducted extensive ablation studies to validate the effectiveness of each component in contributing to the overall performance.", "weaknesses": "- There is a limited set of hand-designed expression forms, which restricts the applicability of the method to other domains.\n- It is a bit illusive how SR can be applied to practical domains.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pnYUHhZNnj", "reviewer": "Reviewer_d71r", "summary": "This paper studies the problem of discovering math equations in real-world complex systems. The authors propose a Reinforcement Symbolic Regression Machine (RSRM) that masters the capability of uncovering complex math equations from scarce data.  Experiments demonstrate the proposed method outperforms baselines on various benchmarks.", "questions": "Please refer to Weaknesses for my questions.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1.\tThe paper is well-written and easy to follow.\n2.\tExperiments demonstrate the proposed method outperforms baselines on various benchmarks.", "weaknesses": "1.\tThe novelty of the proposed method is unclear. The proposed method seems to be a simple combination of existing methods, including the monte carlo tree search (MCTS) algorithm, double q-learning method, and genetic programming method.\n2.\tThe motivation and advantages of the proposed method over previous work are unclear.\n3.\tThe authors claim that they use a MCTS agent for exploration. However, I found that they use greedy selection rather than exploration based on the upper confidence bound, which is a key component for efficient exploration in MCTS. Thus, whether the MCTS agent can explore the environment efficiently is unconvincing. \n4.\tThe proposed modulated sub-tree discovery incorporate three specific search forms into the algorithm. However, the three forms may be too specific to be generally applicable to complex real-world problems.\n5.\tThe authors evaluate the generalization ability of their method in a toy environment. It would be more convincing if the authors could conduct the generalization experiments on more complex real-world problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VwY0oUXLIq", "reviewer": "Reviewer_AFH1", "summary": "This paper proposes a new symbolic regression discovery method that uses Monte Carlo Tree Search for exploration, Double Q-Learning for exploitation, and a sub-tree discovery block to capture equation decompositions. All three components combined in their method named Reinforcement Symbolic Regression Machine, appear to achieve state-of-the-art results on a diverse range of standard symbolic regression benchmarks.", "questions": "* Can you define $c$ in the UCT algorithm in the text when you introduce it?\n* Can you provide empirical evidence for the statement of \"we find that Gaussian random numbers with a unit mean and variance provide more effective initial values for optimization.\", perhaps in an additional appendix?\n* Can you highlight why only 12 AI Feynman equations were used and possibly include results for all the AI Feynman equations in the AI Feynman problem set?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* The proposed approach appears novel and is well-placed in the relevant literature. The introduction does an excellent job of introducing the related work in a concise setting.\n* The paper is well-written and easy to follow.\n* The approach was empirically demonstrated to achieve state-of-the-art performance across the existing SR benchmark tasks and problem sets.\n* The appendix is detailed and extensive, adding further empirical evidence to the core claims.", "weaknesses": "* There are no theoretical results to verify the proposed approach. However, the method does obtain strong empirical performance.\n* Page 3, top. $c$ is undefined. It could perhaps be helpful to state what $c$ is in the UCT algorithm or provide a reference to either another paper or a forward reference to where it is defined.\n* Page 5, \"we find that Gaussian random numbers with a unit mean and variance provide more effective initial values for optimization.\". There seems to be no empirical evidence in the paper supporting this. Although you may have results for this, it could be helpful for the reader if you included an additional Appendix experiment and linked it to this statement as the empirical evidence for this.\n* There are no error bars for the main results tables, Table 1, Table 2, and Table 3.\n* In Table S7, only 12 AI Feynman equations are used, whereas the original AI Feynman paper proposes 100 AI Feynman equations. Why were these 12 AI Feynman equations used, and is it possible to provide results for the complete set of AI Feynman equations?\n\n\nTypos:\n* Page 5: Splitting by Addition paragraph: \"Expand\" -> \"Expanding\"\n* Page 5: Splitting by Addition paragraph: \"Split\" -> \"Splitting\"\n* Page 5: Splitting by Addition paragraph: \"converts\" -> \"convert\"\n* Page 5: Splitting by Addition paragraph: \", then transforms\" -> \", and then transforms\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "dFxZjfh4pp", "reviewer": "Reviewer_XG8a", "summary": "The paper proposes a novel Reinforcement Symbolic Regression Machine (RSRM) method for symbolic regression. The method is based on Monte Carlo tree search (MCTS) and double Q-learning. The method is evaluated on several benchmark datasets and compared with several representative baseline models. The paper contains a detailed ablation study to demonstrate the effectiveness of each component of the proposed method.", "questions": "- (Abstract) \"In nature, the behaviors\". Consider changing to \"In nature, the behavior\"\n- (Abstract) \"Automatically distilling these equations from limited data is cast as a symbolic regression (SR)\". This sentence is confusing. Symbolic Regression **is** the task of distilling equations from data. The sentence gives the impression that SR is something else and that the authors are using SR to solve that something else. Please clarify.\n- (Abstract) \"The discrete search space tends toward infinity\". This statement is not mathematically rigorous. The search space in SR (with variable length expressions) is always infinite. Do the authors mean the complexity of the search space? How do they measure it? Please clarify.\n- (Introduction) \"The early process of fitting expressions involves\". What do the authors mean by \"early process\"? Are they taking a historical perspective of the field? If so, please clarify.\n- (Introduction) \" By incorporating double Q- learning into MCTS, we effectively address issues of overfitting and generate new expressions\". How does Q-learning address overfitting? Please clarify. What do the authors mean with \"generate new expressions\"? Any other SR method is \"generating new expressions\" as well.\n- (Method) \" we introduce an interpolation method (e.g., data pre-processing)\". What do the authors mean by \"interpolation method\" and how is \"data pre-processing\" related to it? Please clarify.\n- (Method) \"the MSDB examines whether the few expressions that perform well adhere to a specific form. For instance, if both $\\exp(x) − x$ and $\\exp(x) + x$ yield favorable results, the expression can be confirmed as $\\exp(x) − f(x)$, thereby allowing us to focus on finding $f(x)$\" How do the authors determine what \"perform well\" means? I can understand the intuition of the authors: if both $\\exp(x) − x$ and $\\exp(x) + x$ yield favorable results, then the expression must be close to the middle of both expressions, i.e., $\\exp(x)$. However, it seems to me that how to determine the threshold in the fitting error to determine what \"perform well\" means is not clear. Is this threshold the same for all datasets? If not how is it determined? Please clarify.\n- (Method) \"This approach effectively reduces the difficulty associated with specific equations.\" This statement is not clear. Please clarify.\n- (Expression Tree) \"The underlying objective of SR is to transform the given task into the generation of an optimal expression tree\". This statement is not clear. Which task is SR transforming into the generation of an optimal expression tree? Please clarify.\n- (Modulated sub-tree discovery) \"This search form focuses on identifying expressions of the form like ex − x and ex + x.\" What do the authors mean by \"search form\"? Please clarify.\n- (Modulated sub-tree discovery) \"In this search form, we obtain good expressions such as 1.57ex and 1.56ex + x,\". The concept of \"search form\" and the use of \"..in this search form\" is not clear. Is search form a synonym of \"search space\"? Please clarify.\n- (Modulared sub-tree discovery) \"The complete from-discovery algorithm, \". There is a typo in \"from-discovery\". Please fix.\n- (Results) How many expressions were evaluated in each algorithm? Is this the same number for all algorithms? This is very important to assess the computational complexity of the algorithm proposed. Please clarify.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The paper contains extensive experiments on a variety of benchmark datasets, demonstrating the superiority of the proposed method over several representative baseline models.\n- The paper contains a detailed ablation study to demonstrate the effectiveness of each component of the proposed method.", "weaknesses": "- The paper has space for improvement in terms of presentation. The paper is not well written and the authors should carefully proofread the paper to improve the presentation quality.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)