Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients
Authors: Brenden K Petersen
Mikel Landajuela Larma
Terrell N. Mundhenk
Claudio Prata Santiago
Soo Kyung Kim
Joanne Taery Kim
Venue: ICLR 2021 Oral (2021)
Keywords: symbolic regression
reinforcement learning
automated machine learning

Abstract: Discovering the underlying mathematical expressions describing a dataset is a core challenge for artificial intelligence. This is the problem of $\textit{symbolic regression}$. Despite recent advances in training neural networks to solve complex tasks, deep learning approaches to symbolic regression are underexplored. We propose a framework that leverages deep learning for symbolic regression via a simple idea: use a large model to search the space of small models. Specifically, we use a recurrent neural network to emit a distribution over tractable mathematical expressions and employ a novel risk-seeking policy gradient to train the network to generate better-fitting expressions. Our algorithm outperforms several baseline methods (including Eureqa, the gold standard for symbolic regression) in its ability to exactly recover symbolic expressions on a series of benchmark problems, both with and without added noise. More broadly, our contributions include a framework that can be applied to optimize hierarchical, variable-length objects under a black-box performance metric, with the ability to incorporate constraints in situ, and a risk-seeking policy gradient formulation that optimizes for best-case performance instead of expected performance.

Published as a conference paper at ICLR 2021
DEEP SYMBOLIC REGRESSION:
RECOVERING MATHEMATICAL EXPRESSIONS FROM DATA VIA RISK-SEEKING POLICY GRADIENTS
Brenden Petersen"
Mikel Landajuela Larma
Lawrence Livermore National Laboratory Livermore, CA, USA bp@llnl.
Lawrence Livermore National Laboratory
Livermore, CA, USA
landa jaatjajeeall@lmme
T. Nathan Mundhenk
Claudio P. Santiago
Lawrence Livermore National Laboratory Livermore, CA, USA
Lawrence Livermore National Laboratory Livermore, CA, USA
mundhenk1@11mg gov
Soo K. Kim
Joanne T. Kim
Lawrence Livermore National Laboratory Livermore, CA, USA
Lawrence Livermore National Laboratory Livermore, CA, USA kim102ellnI. gov
kim7901lnl.
ABSTRACT
Discovering the underlying mathematical expressions describing a dataset a core challenge for artificial intelligence. This is the problem of symbolic regression. Despite recent advances in training neural networks to solve complex tasks, deep learning approaches to symbolic regression framework that leverages deep learning for symbolic regression via a simple idea: large model to search the space of small models. Specifically, use a recurrent neural network to emit a distribution over tractable mathematical expressions and employ novel risk-seeking policy gradient train the network to generate better-fitting expressions. Our algorithm outperforms several baseline methods (including Eureqa, the gold standard for symbolic regression) in its ability to exactly recover symbolic expressions on a series of benchmark problems, both with and without added noise. More broadly, our contributions include framework that can be applied to optimize hierarchical, variable-length objects under black- box performance metric, with the ability to incorporate constraints in situ, and a risk-sceking policy gradient formulation that optimizes for best-case performance instead of expected performance.
underexplored. We propose
INTRODUCTION
Understanding the mathematical relationships among variables in physical system is an integral component the scientific process. Symbolic regression aims to identify these relationships by searching over the space of tractable (i.e. concise, closed-form) mathematical expressions to best fit a dataset. Specifically, given a dataset (X, y). where each point X; E RR" and ys E R, symbolic regression aims to identify function f I R that best fits the dataset, where the functional form of f is a short mathematical expression. The resulting expression can be readily interpreted and/or provide useful scientific insights simply by inspection. In contrast, conventional regression imposes a single model structure that is fixed during training, often chosen to be expressive (e.g. a neural network) at the expense of being easily interpretable.
"Corresponding author.

---

Published as a conference paper at ICLR 2021
Symbolic regression exhibits several unique features that make it an excellent test problem for benchmarking automated machine learning (AutoML) and program symthesis methods: (1) there exist well-established, challenging benchmark problems with stringent success criteria (White et al., 2013);
(2) there exist well-established baseline methods (most notably, the Eureqa algorithm (Schmidt & Lipson )009:; and (3) the reward function is computationally expedient, allowing sufficient experiment replicates to achieve statistical significance. Most other AutoML tasks, e.g. neural architecture search (NAS), do not exhibit these features; in fact, even simply evaluating the efficiency of the discrete search itself is known challenge within NAS (Yu et al. 2019).
The space of mathematical expressions is discrete (in model structure) and continuous (in model parameters), growing exponentially with the length of the expression, rendering symbolic regression a challenging machine learning problem-
to be NP-hard (Lu et al., (016). Given this large, combinatorial search space, traditional approaches to symbolic regression typically utilize evolutionary algorithms, especially genetic programming (GP) (Koza 1992 Schmidt & Lipson 2009). Bäck et al. 2018). In GP-based symbolic regression, population of mathematical expressions is "evolved" using evolutionary operations like selection, crossover, and mutation to improve a fitness function. While GP can be effective, it is also known to scale poorly to larger problems and to exhibit high sensitivity to hyperparameters.
Deep learning has permeated almost all areas of artificial intelligence, from computer vision (Krizhevsky et al.] 2012) to optimal control (Mnih et al.] (015). However, deep learning may seem incongruous with even antithetical toward symbolic regression, given that neural networks are typically highly complex, difficult to interpret, and rely on gradient information. We propose a framework that resolves this incongruity by tying deep learning and symbolic regression together with a simple idea: use a large model (i.e. neural network) to search the space of small models (i.e. symbolic expressions). This framework leverages the representational capacity of neural networks to generate interpretable expressions, while entirely bypassing the need to interpret the network itself.
We present deep symbolic regression (DSR), a gradient- based approach for symbolic regression based on reinforcement learning. In DSR, a recurrent neural network (RNN) emits a distribution over mathematical expressions. Expressions are sampled from the distribution, instantiated, and evaluated based on their fitness to the dataset. This fitness is used as the reward signal to train the RNN using novel risk-seeking policy gradient algorithm. As training proceeds, the RNN adjusts the likelihood of an expression relative to its reward, assigning higher probabilities to better expressions.
We demonstrate that DSR outperforms several baseline methods, including two commercial software algorithms. We summarize our contributions follows: (1) novel method for symbolic regression that outperforms several baselines on a set of benchmark problems, (2) an autoregressive generative modeling framework for optimizing hierarchical, variable-length objects that accommodates in situ constraints, and (3)= a novel risk-seeking policy gradient objective and accompanying Monte Carlo estimation procedure that optimizes for best-case performance instead of average performance.
RELATED WORK
Deep learning for symbolic regression. Several recent approaches leverage deep learning for symbolic regression. AI Feynman (Udrescu & Tegmark| 2020) propose a problen-simplificttion tool for symbolic regression. They use neural networks to identify simplifying properties in dataset (e.g. multiplicative separability, translational symmetry), which they exploit to recursively define simplified sub-problems that can then be tackled using any symbolic regression algorithm. In Grammar VAE, Kusner et al. (2017) develop a generative model for discrete objects that adhere to a pre-specified grammar, then optimize them in latent space. They demonstrate this can be used for symbolic regression; however, the method struggles to exactly recover expressions, and the generated expressions are not always syntactically valid, (Shoo et al. (2018) develop a symbolic regression framework using neural networks whose activation functions are symbolic operators. While this approach enables an end-to-end differenniable system. backpropagation through activation functions like division or logarithm requires the authors to make several simplifications the search space, ultimately precluding learning certain simple classes of expressions like vz or sin (z/v). We address and/or directly compare to these works in Appendices C and
2

---

Published as a conference paper at ICLR 2021
AutoML and program programssytthesis. Symbolic regression is related to both automated machine learning (AutoML) and program synthesis, in that they all involve a search for an executable program (i.e. expression) to solve a particular task (i.e. to fit data) (Abolafia et al. 2018: Devlin et al., [017]] Riedel et al. 2016). More specifically, our framework has many parallels to a body of works within A that use autoregressive RNN define distribution over discrete objects and use reinforcement learning to optimize this distribution under a black-box performance metric (Zoph) & Le 2017: Ramachandran al. 2017 Bello al., 2017, Abolafia et al. 2018). For example, in neural architecture search (Zoph & Le 2017). , an RNN searches the space of neural network architectures, encoded by a sequence of discrete "tokens" specifying architectural properties (e.g. number of neurons) of each layer. The length of the sequence is fixed or scheduled during training: in contrast, our framework defines a search space that both inherently hierarchical and variable length. Ramachandran et al. (2017) search the space of neural network activation functions. While this space is hierarchical in nature, the authors (rightfully) constrain it substantially by positing a functional unit that is repeated sequentially. thus restricting their search space back to fixed-length sequence. However, a repeating-unie constraint is not practical for symbolic regression because the ground truth expression may have arbitrary structure.
Autoregressive models. The Raseb distribution over expressions used in DSR is autoregressive, meaning each token is conditioned on the previously sampled tokens. Autoregressive models have proven to be useful for audio and image data (Oord et al.] 2016a bb) in addition to the AutoMI. works discussed above; we further demonstrate their efficacy for hierarchical expressions. GraphRNN defines a distribution over graphs that generates an adjacency matrix one column at a time in autoregressive fashion YYou al. 2018). In principle, GraphRNN could be constrained to define distribution over expressions, since trees are a special case of graphs. However, GraphRNN constructs graphs breadth-firss whereas expressions are more naturally represented using depth- first traversals (Li et al. (005). Further, DSR exploits the hierarchical nature of trees by providing the parent and sibling as inputs to the RNN, and leverages the additional structure of expression trees that node's value determines its number of children (e.g. cosine is unary operator and thus has one child).
Risk-aware reinforcement learning. Many of the AutoML methods discussed above suffer from what we call the "expectation problem." That is, policy gradient methods are fundamentally suited for optimizing expectations; however, domains like neural architecture search and symbolic regression evaluated by the few single best-perfor rest-performing samples. Thus, there is a mismatch between the training objective function and the true desired objective: to maximize best-cas best-case performance. Abolatia al. (2018) address the expectation problem by maintaining a priority queue of the best seen samples and using supervised learning to increase the likelihood of those top samples. Similarly, Liang et al. (2018) use a memory buffer to augment a policy gradient with off off-policy training. These methods, however, only apply in the context of reinforcement learning environments with both deterministic transition dynamics and deterministic rewards. In contrast, the risk-seeking policy gradient introduced here is general, applying to any reinforcement learning environment and any stochastic policy gradient algorithm trained using batches, e.g. on Atari using proximal policy optimization (Schulman et al., 2017]. Lastly, our risk-seeking policy gradient is closely related to the EPOpt-S algorithm used for robust reinforcement learning (Rajeswaran et al. 2016), which is based on risk-avers, policy gradient formulation (Tamar et al. 2014).
METHODS
Our overall approach involves representing mathematical expressions as sequences, developing an autoregressive model to generate expressions under pre-specified set of constraints, and developing a risk-seeking policy gradient to train the model to generate better-fitting expressions.
3.1 GENERATING EXPRESSIONS WITH A RECURRENT NEURAL NETWORK
We leverage the fact that mathematical expressions can be represented using symbolic expression trees. Expression trees are a type of binary tree in which internal nodes are mathematical operators and terminal nodes are input variables or constants. Operators may be unary (i.e. one argument, such as sine) or binary (i.e. two arguments, such as multiply). Further, as a sequence of node values or "tokens" by using its pre-order traversal (i.e. by visiting each node depth-first, then left-to-righh). This allows us to generate an expression tree sequentially while still
represent an expression tree
3

---

Published as a conference paper at ICLR 2021
Pre-order traversal
Library
Sampled
token:
Categorical distribution:
Expression free
RNN:
sin
log
Parent:
Sibling:
Figure 1: A. Example of sampling expression from the RNN. For each token, the RNN emits categorical distribution over tokens, a token sampled, and the parent and sibling of the next token used as the next input to the RNN. Subsequent tokens sampled autoregressively until the tree is complete (i.e. all tree branches reach terminal nodes). The resulting sequence tokens is the tree's pre-order traversal, which can be used reconstruct the and instantiate its corresponding expression. Colors correspond to the number of children for each token. White circles represent empty tokens. Numbers indicate the order in which tokens were sampled. B. The library of tokens. C. The expression tree sampled in A. In this example, the sampled expression is sin(cr), (og(y), where the value of the constant is optimized with respect to an input dataset.
maintaining one-to-one correspondence between the tree and its traversal Thus, we represent an expression by the pre-order traversal of its corresponding expression tree We denote the 1th token of the traversal as T; and the length of the traversal as |T| = T. Each token has a value within a given library L of possible tokens, e.g. {+,
+, s,,in, x/.
We generate expressions one token at time along the pre-order traversal (from T1 to TT). categorical distribution with parameters & defines the probabilities of selecting each token from L. To capture the "context" of the expression as is being generated, we condition this probability upon the selections of all previous tokens in that traversal. This conditional dependence can be achieved very generally using an RNN with parameters 0 that emits probability vector zó in an autoregressive manner. Specifically, the 7th output of the RNN passes through softmax layer (with shared weights across time steps) produce vector 26(*), which defines the probability distribution for selecting the (4)
7th token r;, conditioned on the previously selected tokens T1:(i-1)- That is, p(T:///: 0) vi D(r.)
where C(T;) is the index in C corresponding to The likelihood of the entire sampled expression is
simply the product of the likelihoods of its tokens: p(T [[]=] p( 8) I [r vi ((T.)
An example of the sampling process is illustrated in Figure [ pseudocode provided in Algorithm in Appendix A Note that different samples from the distribution have different tree structures of different size; thus, the search space is inherently both hierarchical and variable length.
Providing hierarchical inputs to the RNN. Conventionally, the input to the RNN when sampling a token would be a representation of the previously sampled token. However, the search space for symbolic regression is inherently hierarchical, and the previously sampled token may actually be very distant from the next token to be sampled the expression tree. For example, the fifth and
In general, a pre-order traversal insufficient uniquely reconstruct the tree. However, we know how many children each node has based on its value, e.g. "multiply" is binary operator and thus has two children. pre-order traversal plus the corresponding number of children for each node is sufficient to uniquely reconstruct the tree.
"Given an expression tree. the corresponding mathematical expression is unique: however, given an expres- sion, its expression tree is not unique. For example, x² and x x are equivalent expressions but yield different trees. For simplicity, use somewhat ahusively to refer to expression where it technically refers an expression tree (or equivalently, pre-order traversal).
4

---

Published as a conference paper at ICLR 2021
sixth tokens sampled in Figure 1 are adjacent nodes in the traversal but are four edges apart in the expression tree. To better capture hierarchical information, we provide as inputs to the RNN a representation of the parent and sibling nodes of the token being sampled. We introduce an empty token for cases in which a node does not have a parent or sibling. Pseudocode for identifying the parent and sibling nodes given partial traversal is provided in Subroutine i Appendix A
Constraining the search space. Under our framework, it is straightforward to apply a priori constraints reduce the search space. To demonstrate, we impose several simple, domain-agnostic constraints: (1) Expressions are limited to a pre-specified minimum and maximum length. We selected minimum length of to prevent trivial expressions and maximum length of 30 to ensure expressions are tractable. (2) The children of an operator should not all be constants, as the result would simply be different constant. (3) The child of unary operator should not the inverse of that operator, e.g. log( (x))(())) is not allowed. (4) Descendants of trigonometric operators should not be trigonometric operators, e-g. sin(x + cos(x)) is not allowed because cosine is descendant of sine. While still semantically meaningful, such composed trigonometric operators do not appear in virtually any scientific domain.
We apply these constraints in situ (i.e. concurrently with autoregressive sampling) by zeroing out the probabilities of selecting tokens that would violate a constraint. Pseudocode for this process is provided in Subroutine in Appendix A This process ensures that samples always adhere to all constraints, without rejecting samples post hoc. In contrast, imposing constraints with GP requires rejecting evolutionary operations post hoc (Fortin al., 2012), which be problematic (Craenen it al. 2001), and as we show in our experiments, can actually reduce performance.
Reward function. Once pre-order traversal sampled, instantiate the corresponding symbolic expression and evaluate it with reward function. A standard fitness measure in GP-based symbolic regression is normalized root-meam mean--suuare error (NRMSE), the root-mean-squee -uuare error normalized by the standard deviation of the target values, Ty. That is, given a dataset (X, y) of size n and candidate
expression f, NRMSE 1 V A Ei= f(X.))² To bound the reward function, we apply By squashing function: R(T) = 1/(1 + NRMSE).
Constant optimization. If the library L. includes the constant token, sampled expressions may in- clude several constant placeholders. These can be viewed as parameters of the symbolic expression, which we optimize by maximizing the reward function: e* = arg maxę R[(;; E). using nonlinear optimization algorithm, e.g. BFGS (Fletcher 2013). We perform this inner optimization loop for each sampled expression as part of the reward computation before performing each training step.
3.2 TRAINING THE RNN USING POLICY GRADIENTS
The reward function is not differentiable with respect to 0: thus, we turn to reinforcement learning to train the RNN to produce fetter-fitting expressions. In this view, the distribution over mathematical expressions p(T(0) is like policy, sampled tokens are like actions, the parent and sibling inputs are like observations, sequences corresponding to expressions are like episodes, and the reward is a terminal, undiscounted reward only computed when an expression completes.
Standard policy gradient. We first consider the standard policy gradient objective to maximize (0), defined as the expectation of reward function R(T) under expressions from the policy:
(0) (E>=p(T(0) [A(r)). The standard REINFORCE policy gradient (Williams) 1992) can be used to maximize this expectation via gradient ascent:
V Jadd ( = V E. [R(T)]
= (prr(()) [R( p(]()]
This result allows one to estimate the expectation using samples from the distribution. Specifically,
an unbiased estimate of V aaaa can be obtained by computing the sample mean over a batch of N
sampled expressions T {r (6) L.V
VoJaa
N Σ
This is an unbiased gradient estimate, but in practice it has high variance. To reduce variance, is common subtract a baseline function from the reward. As long as the baseline is not a function of
5

---

Published as a conference paper at ICLR 2021
the current batch of expressions, the gradient estimate is still unbiased. Common choices of baseline functions are a moving average of rewards or an estimate of the value function. Intuitively. the gradient step increases the likelihood of expressions above the baseline and decreases the likelihood of expressions below.
Risk-seeking policy gradient. The standard policy gradient objective, (0), is defined as an expectation. This is the desired objective for control problems in which one seeks to optimize the average performance of policy. However, in domains like symbolic regression, program synthesis, or neural architecture search, the final performance is measured by the single or few best-perforning samples found during training. Similarly, one might be interested in a policy that achieves "high score" in control environment (e.g. Atari). For such problems, (s) is not an appropriate objective, as there is mismatch between the objective being optimized and the final performance metric: this is the "expectation problem." 'oo address this disconnect, we propose an alternative objective that focuses learning only on maximizing best-case performance. We first define Rg as the (1 e--quntiele of the distribution of rewards under the current policy. Note this is function of 0 but is typically intractable. We then propose a new learning objective, Jria (0; e), parameterized by e:
(0; = [T-p((10))
[R(T) (((T) 2 Rg (0)]
(1)
This objective aims increase the reward of the top fraction of samples from the distribution, without regard for samples below that threshold. This objective bears close resemblance with e- conditional value at risk (CVaR), for which the "s" symbol "ymbol is used instead of "2" and the E-quantile of rewards is used instead of the (1 - e'-quantile. Optimizing CVaR is a form of risk-averse learning that results in policy that is robust against catastrophic outcomes (Tamar al. 2014 Rajeswaran et al. 2016). In contrast, optimizing (nss(0;;) yields a risk-seeking policy gradient that aims to increase best-case performance the expense of lower worst-case and average performances. Next, we show the analogous policy gradient of Jaix (0;=) and how to estimate it via Monte Carlo sampling. Proposition 1. Let Jri (0;e) denote the conditional expectation of rewards above the (1 (1-8) (1-s)-quannie -quantile, as in Equation 7. Then the gradient of Jrisk =):=) is given by:
V o Jris (0;=) = Er-p(r (r)))(((0)))))))))))))))) (0) [R(T) - (0)) V log R(T) 2 (0)
The proof and assumptions are provided in Appendix and are adapted from the policy gradient
derivation for the CVaR objective (Tamar et al. 2014). This result suggests simple Monte Carlo
estimate of the gradient from batch of samples:
VaJra (0;E) R
E [(( - R=((0] n(- (0) V log logp(>(")0),
where Rg (0) is the empirical (1 11-53 =) -uquntile of the batch of rewards, and 1g returns if condition is true and 0 otherwise. Essentially, this is the standard REINFORCE Monte Carlo estimate with two differences: (1) theory suggests specific baseline, R,((),, whereas the baseline for standard policy gradients is non-specific, chosen by the user; and (2) effectively, only the top E fraction of samples from each batch are used in the gradient computation. This process is essentially the opposite of
the approach used to optimize CVaR (Tamar et al. 2014 Rajeswaran et al. 2016) for risk-averse reinforcement learning, in which only the bottom ε fraction of samples from each batch are used.
Note that in Proposition| I and the corresponding Monte Carlo estimation procedure, 7 need not refer symbolic expression or even discrete object. For example, it may refer to a state-action trajectory in a control problem. Thus, the risk-secking policy gradient formulation is general and can easily be applied to any reinforcement learning environment using any stochastic policy gradient algorithm trained using batches, e.g. proximal policy optimization (Schulman et al.) (017).
Lastly, in accordance with the maximum entropy reinforcement learning framework (Haarnoja et al. 2018), we provide a bonus to the loss function proportional to the entropy of the sampled expressions. Pseudocode for DSR is shown in Algorithm [ Source code is made available at https://githuu. htw.aa/.cconooeeeennttteeeeeeeee
deep- regresssion.
4 RESULTS AND DISCUSSION
Evaluating DSR. We evaluated DSR on the Nguyen symbolic regression benchmark suite (Uy et al. 201 set of 12 commonly used benchmark expressions developed and vetted by the symbolic

---

Published as a conference paper at ICLR 2021
Algorithm 1 Deep symbolic regression with risk-seeking policy gradient
input learning rate at entropy coefficient Aæς; risk factor E; batch size N; reward function R output Best fitting expression 7*
I: Initialize RNN with parameters 8, defining distribution over expressions p(-10)
2: repeat
3: T f- {-(0) p(-10));=1
Sample N expressions (Alg, 2 in Appendix Ap Optimize constants w.r.t. reward function
4: T - (OpiiiizeConstattsrr))) R))
5:
R < {
> Compute rewards D. Compute reward threshold > Select subset of expressions above threshold
6:
RS < (1 (1-6 8) e--unnielle of R
7:
T t (r()) R(+(0) 2 R
8:
R + (R(+()): R()()) 2 R
b Select corresponding subset of rewards Compute risk-seeking policy gradient D Compute entropy gradient > Apply gradients
9:
91 ReduceMean (R -R,) Vo (ogg(((0))
10:
92 + ReduceMean Ay (oH(T
0 0 + a(gi 92)
if maax R [((T)) then T* (argmax max R)
b Update best expression
13: return T*
regression community (White et al. 2013). Each benchmark is defined by a ground truth expression, training and test dataset, and set of allowed operators, described in Table in Appendix [D] The training data is used to compute the reward for each candidate expression, the test data is used to evaluate the best found candidate expression the end of training, and the ground truth expression is used to determine whether the best found candidate expression was correctly recovered.
We compared DSR against five strong symbolic regression baselines: (1) PQT:
an implementation of our framework trained using priority queue training (Abolatia et al. (0188) in place of the risk-seeking policy gradient, (2) VPG: a "vanilla" implementation of our framework using the standard policy gradient (with baseline) in place of the risk-seeking policy gradient, (3) GP: a standard GP-based symbolic regression implementation (Fortin et al. 2012). (4) Eureqa: popular commercial softward based on (Schmidt& & Lipson (2009), the gold standard for symbolic regression, and (5) Wolfram: commercial baased on Markov chain Monte Carlo and nonlinear regression (Fortuna 2015). Each baseline is further detailed in Appendix C Notably, the two RNN-based baselines (PQT I' and VPG) differ from DSR only via their training objective, i.e. a supervised objective for PQT and the standard policy gradient for VPG. All other aspects of our framework (in situ constraints, hierarchical RNN inputs, constant optimization) are also included for these baselines. While GP can incorporate post hoc constraints, we found that doing so actually hinders performance (see Appendix [F); thus, they are excluded for GP we exclude all constraints other than the maximum length constraint.
Each experiment consists of million expression evaluations for DSR, PQT, VPG, and GP, by which point training curves have levelled off. Hyperparameters were tuned on benchmarks Nguyen-7 and Nguyen-10 using grid search comprising 800 hyperparameter combinations for GP and 81 combinations for each of DSR, PQT, and VPG (see Appendix D for details). Commercial software algorithms (Eureqa and Wolfram) do not expose hyperparameters and were run until completion. All experiments were replicated with 100 different random seeds for each benchmark expression. Wolfram results are only presented for one-dimensional benchmarks because the method is not applicable higher dimensions. Additional experiment details provided in Appendix D Training curves are provided in Appendix F
In Table I report the recovery rate for each benchmark. We use the strictest definition of recovery: exact symbolic equivalence, as determined using a computer algebra system, e.g. SymPy (Meurer al. 2017). In Table Appendix F report recovery on several additional variants of Nguyen benchmarks in which we introduced real-valued constants (to demonstrate the constant optimizer) or altered the functional form to make the problems more challenging. DSR significantly outperforms
all five baselines in its ability to exactly recover benchmark expressions. In Tables 5 8 in Appendix
E we compare DSR to literature-reported values from four additional baseline methods (Huynh)
et al. 2016
Kusner et al. 2017 Jin et al. 2019 Trujillo et al. 2016) by carefully recapitulating their
*Sold by DataRobot, Inc. (www.datarobot.com comb and formerly by Nutonian, Inc. (www.nutonian.com) 'Sold by Wolfram Research, Inc. as part of Mathematica www.wolfram.ccom/aaaeaaiiaio
7

---

Published as a conference paper at ICLR 2021
Table 1: Recovery rate comparison of DSR and five baselines on the Nguyen symbolic regression benchmark suite. A bold value represents statistical significance (p 10 3) across all benchmarks.
Benchmark
Expression
DSR
PQT
VPG
GP
Eureqa
Wolfram
N
x
+ x" +x
100%
100%
96%
100%
100%
100%
Nguyen-2
x t x + x² to
100%
99%
47%
97%
100%
100%
N
+ x x te² 2 tx x²
100%
86%
4%
100%
95%
100%
Nguyen-4
x + x² +
100%
93%
100%
70%
100%
Nguyen-5
sin(x ( coo((2) -1
72%
73%
5%
45%
73%
2%
N
sin(r) sin(x + + log(æ + 1) + log(a" +1
100%
98%
100%
91% %
100%
1%
Nguyen-7
35%
41%
3%
0%
85%
0%
Nguyen-8
V sin(x) sin(y²)
96%
21%
5%
5%
0%
71%
Nguyen-9
100%
100%
100%
100%
100%
Nguyen-10 10
2 2sin(x) cos(y) x²
100%
91%
99%
76%
64%
Ng Nguyen-11
100%
100%
100%
7%
100%
Nguyen-12
x²
+ 20
y
0%
0%
0%
0%
0%
Average
83.6%
75.2%
46.7%
60.1%
73.9%
10
Risk-soeking
Risk-socking Top batch
Densty
Full batch
0.6
0
2k 0.4
10
Standard PG, Full batch
D
Standard PG, Top batch
Densky
Full
10
Mide saokingiin PG
Top
Best
0.0
Reward
Reward
Expressions
Figure 2: A - D. Empirical reward distributions for Nguyen-8. Each curve is Gaussian kernel density estimate (bandwidth 0.25) of the rewards for a particular training iteration, using either the full batch of expressions (A and C) or the top e fraction of the batch (B and D), averaged over all training runs. Black plots and B) were trained using the risk-seeking policy gradient objective. Blue plots (C and D) were trained using the standard policy gradient objective. Colorbars indicate training step. Triangle markings denote the empirical mean of the distribution the final training step. E. Training curves for mean reward of full batch (dotted), mean reward of top fraction of the batch (dashed), and best expression found so far (solid), averaged over training runs.
experimental setup (e.g. choice of library, benchmarks, and measure of performance). DSR greatly outperforms each study's published results.
Characterizing the risk-seeking policy gradient. The intuition behind the risk-sccking policy gradient is that it explicitly optimizes for best-case performance, possibly at the expense of average performance. We demonstrate this visually in Figure by comparing the empirical reward distribu- tions when trained with either the risk-seeking or standard policy gradient for Nguyen-8. (Analogous plots for all Nguyen benchmarks are provided in Appendix F Interestingly, at the end of training, the mean reward over the full batch (an estimate of Jau (0) larger when training with the standard policy gradient, even though the risk-seeking policy gradient produces larger mean over the top fraction of the batch (an estimate of Jrsk(B (6; E)) and superior best expression. This is consistent with the intuition of maximizing best-case performance at the expense of average performance. In contrast, the best-case performance of the standard policy gradient plateaus early in training (Figure 2E, dashed blue curve), whereas the risk-secking policy gradient continues to increase until the end of training (Figure 2E, dashed black curve).
Ablation studies. Algorithm includes several additional components relative to "vanilla" policy gradient search. We performed a series of ablation studies to quantify the effect of each of these components, along with the effects of the various constraints the search space. In Figure 3 show recovery rate for DSR on the Nguyen benchmarks for each ablation. While no single ablation leads to catastrophic failure, combinations of ablations can cause large degradation in performance.
8

---

Published as a conference paper at ICLR 2021
No ablations
80%
data
10x data
100%
No entropy constraints No risk/ent. risk/con.
Recovery
600%
Recovery
75%
40%
50%
No ent./con.
+ DSR + PQT
+ VPG
25%
DSR
ablations
GP
Eurega
PQT
0%
40%
80%
0.02
0.04
0.06
0.08
0.1
0.0
0.1 0.2 0.3 0.4 0.5 Reward noise
Recovery
Noise level
Figure 3: Recovery for various ablations of Algorithm across all Nguyen benchmarks. Error bars represent standard error.
Figure Recovery dataset noise and dataset size across all Nguyen benchmarks. Error bars represent stan- dard error.
Figure 5: Recovery vs added reward noise on Nguyen-4. Error bars rep- resent standard error.
Noisy data and amount of data. We evaluated the robustness of DSR to noisy data by adding independent Gaussian noise the dependent variable, with mean zero and standard deviation proportional to the mous-mean-suaaee of the dependent variable the training data. In Figure 4 varied the proportionality constant from 0 (noiseless) to 10-1 and evaluated each algorithm (except Wolfram, which catastrophically fails for even the smallest noise level) across all Nguyen benchmarks. Because expressions can overfit the noise in these experiments, as exact symbolic equivalence on any expression along the reward cwwrl--omplexxtty Pareto front at the end of training (see Appendix Difor details). Increasing the dataset size may help prevent overfitting by smoothing the reward function. Thus, we repeated the noise experiments with 10-fold larger training datasets (Figure 4 dashed lines). For each dataset size, DSR consistently outperforms all baselines.
defined recovery
Performance under reward noise. The resk-ngeking policy gradient is derived from a conditional expectation, and is thus well-justified for tasks with stochastic rewards. In contrast, PQT assumes deterministic rewards. Since symbolic regression is deterministic task, we emulated a stochastic reward function by adding independent Gaussian noise directly to the reward function: R'(r) = R(T) + (0, o). In Figure [ we compare DSR and PQT performance under increasing reward noise on benchmark Nguyen-4, a task that an easier exploitation problem but difficult exploration problem. As reward noise increases, recovery for PQT heavily relies exploration and luck, whereas DSR recovery remains stable even for high reward noise.
5 CONCLUSION AND FUTURE WORK
introduce a reinforcement learning approach to symbolic regression that outperforms state-of-the- art baselines in its ability recover exact expressions benchmark tasks. Since both DSR and GP generate expression trees, there are many opportunities for hybrid methods, for example including several generations of evolutionary operations within the inner optimization loop. Our framework includes a flexible distribution over hierarchical, variahle-length objects that allows imposing in situ constraints. Our framework is easily extensible other domains, which save for future work; for example, searching the space of expressions to be used as control policies in reinforcement learning environments, searching the space of organic molecular structures for high binding affinity to reference compound. Our risk-seeking policy gradient formulation can also be applied to more traditional reinforcement learning domains; for example, optimizing for high score (instead of average score) in Atari video games.
ACKNOWLEDGMENTS
We thank Ruben Glatt, Thomas Desautels, Priyadip Ray, David Widemann, and the 2019 UC Merced Data Science Challenge participants for their useful comments and insights. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NNA2344 Lawrence Livermore National Security, LLC. LLNI-CCONF-700457.

---

Published as a conference paper at ICLR 2021
REFERENCES
Daniel A Abolatia, Mohammad Norouzi, Jonathan Shen, Rui Zhao, and Quoc V Le. Neural program synthesis with priority queue training, arXiv preprint arXiv: 00056,, 2018.
Thomas Bäck, David B Fogel, and Zbigniew Michalewicz. Evolutionary Computation 1: Basic Algorithms and Operators. CRC press, 2018.
Irwan Bello, Barret Zoph, Vijay Vasudevan, and Quoe V Le. Neural optimizer search with reinforce- ment learning. In Proceedings the 34th International Conference on Machine Learning- -oolue 70, pp. 459-468. JMLR. org, 2017.
BGW Craenen, AE Eiben, and Marchiori. How to handle constraints with evolutionary algorithms. Practical Handbook Of Genetic Algorithms: Applications, pp. 341-361, 2001.
Jacob Devlin, Jonathan Uesato, Surya Bhupatiraju, Rishabh Singh, Abdel-rahman Mohamed, and Pushmeet Kohli. Robustfill: Neural program learning under noisy ilo. arXiv preprint arXiv:100..07469, 2017.
Roger Fletcher. Practical Methods of Optimization. John Wiley & Sons, 2013.
Félix-Antoine Fortin, François-Michel De Rainville, Marc-André Gardner, Marc Parizeau, and Christian Gagné. Deap: Evolutionary algorithms made easy. Journal of Machine Learning Research, 130ul):2111-2175 2012.
Giorgia Fortuna. Automatic formula discovery in the wolfram language. In Wolfram Technol- ogy Conference 2015. Wolfram Research, 2015. URL hRLhhtts:///:bbrrr wol www.fam........mmmm com/
Infocenterr/Coneereeees 9
Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine. Soft actor-critic: Off-policy maxi- mum entropy deep reinforcement learning with a stochastic actor. arXiv preprint arXiv:1801.01290 2018.
Quang Nhat Huynh, Hemant Kumar Singh, and Tapabrata Ray. Improving symbolic regression through a semantics-driven framework. In 2016 IEEE Symposium Series on Computational Intelligence (SSCI), pp. 1-8. IEEE, 2016.
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression. arXiv preprint 1xXi:19100088/)) 2019.
John R Koza. Genetic Programming: On the Programming of Computers by Means of Natural Selection, volume MIT press, 1992.
Alex Krizhevsky, Iyya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolu- tional neural networks. In Advances in neural information processing systems, pp. 1097-1105, 2012.
Matt J Kusner, Brooks Paige, and José Miguel Hernández-Lobato. Grammar variational autoencoder. In Proceedings of the 34th International Conference on Machine Learning- Volume 70, pp. 1945- 1954. JMLR. org, 2017.
Xin Li, Chi Zhou, Weimin Xiao, and Peter C. Nelson. Prefix gene expression programming. In Late breaking paper at Genetic and Evolutionary Computation Conference (GECCO 2005), Washington, DC, USA, pp. 25-29, 2005.
Chen Liang, Mohammad Norouzi, Jonathan Berant, Quoe Le, and Ni Lao. Memory augmented policy optimization for program synthesis with generalization. arXiv preprint arXiv.1807.02322. 2018.
Qiang Lu, Jun Ren, and Zhiguang Wang. Using genetic programming with prior formula knowledge to solve symbolic regression problem. Computational intelligence and neuroscience, 2016, 2016.
10

---

Published as a conference paper at ICLR 2021
Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondiej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore. Sartaj Singh, Thilina Rathmayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Stěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. Sympy: symbolic computing in python. PeerJ Computer Science, 3:3103, January 2017. ISSN 2376-5992. doi: 10.7717/peerj-os. 103. URL hwt.tp::://o..oooo/oee//o
. 7717/ -cc.:peeeerecc.c.000 103
Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei Rusu, Joel Veness, Mare G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. Human-level control through deep reinforcement learning. Nature, 518(7540):529, 2015.
Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu. Wavenet: A generative model for raw audio. arXiv preprint arXiv:1609.034499 2016a.
Aaron van den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu. Pixel recurrent neural networks. arXiv preprint arXiv: 00069.0...69..... 2016b.
Aravind Rajeswaran, Sarvjeet Ghotra, Balaraman Ravindran, and Sergey Levine. Epopt: Learning robust neural network policies using model ensembles. arXiv preprint arXiv: 2016.
Prajit Ramachandran, Barret Zoph, and Quoe V Le. Searching for activation functions. arXiv preprint arXiv:1710.05941, 2017.
Sebastian Riedel. Matko Bosnjak, and Tim Rocktäschel. Programming with a differentiable forth interpreter. CoRR, abs/1605.06640 2016.
Subham S Sahoo, Christoph H Lampert, and Georg Martius. Learning equations for extrapolation and control. arXiv preprint
2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data, science, 324(5923):81--8) 2009.
John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint 1777v:1707.06347, 2017.
Aviy Tamar, Yonatan Glassner, and Shie Mannor. Policy gradients beyond expectations: Conditional value-at-risk. arXiv preprint arXiv: 44.4.34444.... 2014.
Leonardo Trujillo, Luis Muñoz, Edgar Galván-López, and Sara Silva. neat genetic programming: Controlling bloat naturally. Information Sciences, 333:21-43, 2016.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: physics-inspired method for symbolic regression. Science Advances, 6(16):eaa26311, 2020.
Nguyen Quang Uy. Nguyen Xuan Hoai, Michael O'Neill, Robert McKay, and Edgar Galván- López. Semantically-based crossover in genetic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Machines, 12(2):91-119, 2011. 2011.
David R White, James Medermott, Mauro Castelli, Luca Manzoni, Brian W Goldman, Gabriel Kronberger, Wojciech Jaškowski, Una-May O'Reilly, and Sean Luke. Better gp benchmarks: community survey results and proposals. Genetic Programming and Evolvable Machines, 14(1): 3-29, 2013.
Ronald J Williams. Simple statistical sendient-folloiing algorithms for connectionist reinforcement learning. Machine learning, 8(3-4):229 256, 1992.
Jiaxuan You, Rex Ying, Xiang Ren, William L Hamilton, and Jure Leskovec. Graphmm: Generating realistic graphs with deep auto-regressive models. arXiv preprint arXiv: 1xXiv:/8000 1802.08775 2018.
Kaicheng Yu, Christian Sciuto, Martin Jaggi, Claudiu Musat, and Mathieu Salzmann. Evaluating the search phase of neural architecture search. arXiv preprint arXiv: /90..089000812 2019
Barret Zoph and Quoe V Le. Neural architecture search with reinforcement learning. International Conference Learning Representations, 2017.
11

---

Published as a conference paper at ICLR 2021
APPENDIX A PSEUDOCODE FOR ADDITIONAL ALGORITHMS AND SUBROUTINES
Pseudocode for sampling expression from the RNN. The sampling process in DSR (line 4 of Algorithm 1) is more complicated than typical autoregressive sampling procedures due to applying constraints in situ and providing hierarchical information to the RNN. Thus, provide pseudocode for this process in Algorithmp 2 Within this algorithm, the function Arity(T:) simply returns the arity (number of arguments) of token 7r i.c. two for binary operators, one for unary operators, or zero for input variables or constants,
Algorithm 2 Sampling expression from the RNN
imput RNN with parameters 0: library of tokens L output Pre-order traversal T of an expression sampled from the RNN
1: +
Initialize empty traversal b> Initialize counter for number of unselected nodes Initial RNN input empty parent and sibling b Initialize RNN cell state to zero
2: counter t
3: x =- +-empty|leplly
4: cp 1-
5: for = 1, do
6:
(2b(), ,c,) + (NN((2,c-1;0))
Emit probabilities; update state P. Adjust probabilities
ei(i) + ApplyConstrrinnss 7)
7:
8:
T + Categorical
D. Sample next token Append token to traversal Update number of unselected nodes If expression is complete, return it D Compute next parent and sibling
9:
+
10:
counter + counter + Arity(n.)
11: if counter then return
12:
x t ParentSibling(7)
Additional subroutines. DSR includes several subroutines used when sampling expression from the RNN and during training. In Subroutine we describe the function ParentSibling(1) used in Algorithm 2 which computes the parent and sibling of the next token be sampled. This subroutine uses the following logic. If the final node in the partial traversal is a unary or binary operator, then that node the parent and there is no sibling. Otherwise, the subroutine iterates backward through the traversal until finding node with an unselected child node. That node is the parent and the subsequent node is the sibling.
Subroutine 1 Computing parent and sibling inputs to the RNN
input Partially sampled traversal т
output Concatenated parent and sibling tokens of the next token to be sampled > LLnngth of partial traversal
1: T T'-Ir |r|
2: counter <-
Initialize counter for number unselected nodes
3: if Arity((T)) 0 then
4: parent + Tr
5: sibling empty
6: return sibbling
7: for T, ..,.., do
Iterate backward Update number of unselected nodes
8: counter + counter + Arity(T:) -1
9:
if counter then
10: parent +
sibling f T++1
12:
return parent |sibling
In Subroutine 2 we describe the function ApplyConstraints((p..c C, used in Algorithm 2 which zeros out the probabilities of tokens that would violate any given constraints. Within this subroutine, the user-specific function ViolatesConstrain((? L,) returns TRUE if adding the 7th token of the library L to the partial traversal would violate any user-specified constraints, and FALSE otherwise.
In Subroutine 3] we describe the function OptimizeConstantsst: R) used Algorithm [ which optimizes the placeholder constants & of expression with respect the reward function using
12

---

Published as a conference paper at ICLR 2021
Subroutine 2 Applying generic constraints in situ when sampling from the RNN
input Categorical probabilities th; corresponding library of tokens L; partially sampled traversal output Adjusted categorical probabilities
1: L
Length of library
2: for = :rri=l..,, do
3: if' V'olatesConnsrainn((.... then dh. f 0 If the token would violate constraint, set its
probability to 0
4: x/ f- E.
> Normalize probability vector back to 1
5: return
black-box optimizer, e-g. BFGS. In Algorithm[] T this subroutine is abstracted away into the computa- tion of the reward function; here, we explicitly consider the reward function R(T;E) as parameterized by the placeholder constants. Within this subroutine, the function ReplaceConstants(:: (T) replaces the placeholder constants in the expression with the optimized constants 5*
Subroutine Optimizing the constants of an expression (inner optimization loop)
input Expression with placeholder constants Ç: reward function Rt output Expression r* with optimized constants {* {)::
E*
>> Maximize reward w.r.. constants, e.g. with BFGS D Replace placeholder constants
2: r² x- Replace(Constant((7..8)
3: return 7*
APPENDIX B (PROPOSITION 1)
PROOF OF POLICY GRADIENT FOR RISK RISK-SEEKING OBJECTIVE
Tamar al. (2014) proved the following result for the policy gradient of the e-conditional value
risk objective, (0):
VoJevaR(O) Voll [R((T) R(T) < (o]
= Er>p(tt0) [(R(T) (O)) V logp(T|0) R(T) S (&)
For completeness, we adapt this proof from the CVaR objective to our risk risk-seeking objective given in (1). We emphasize that the proof closely follows [amar et al. (2014). The difference amounts to defining the threshold R. as the top quantile of rewards (instead of the bottom & quantile), and replacing R(T) < R. in the CVaR objective with R(T) 2 R. in our objective. This difference results in the limits of integration swapping, causing additional minus sign that eventually cancels out. Important differences from the proof in Tamar et al. (2014) are colored in red. As in Tamar et al. (2014), first demonstrate the policy gradient in the single variable case, then extend the multivariable case in the reinforcement learning setting. The assumptions follow Assumptions detailed in Tamar et al. (2014).
Proof. Consider a bounded random variable Z E [-b, b] generated from parameterized distribution 0)) The (1 e) quantile of is:
Qz=s(Z;0) 0) = inf(= : CDF(>) 2 1-
where CDF(2) is the cumulative distribution function corresponding to p(Z We define the risk-seeking objective Visk(0;s) as the expectation of the fraction of the best outcomes of Z: (fi E) = Σz-m<z1)
Z 2 - (Z;
We define D3 as the set of values above this quantile:
D, = f2 [ b, b] : 2 (Z:0)
By construction, Dg is simply the interval (Z:0) b). and
[ De
plz10)dz =
(3)
13

---

Published as a conference paper at ICLR 2021
Rewriting the conditional expectation in Equation ( an integral,
Jiik (0: (0:
p(z)0)adz
JeEn, p(210)dz =DDO
1
c p(=10)zdz
8 $JEen
1
p(=10)=dz
J.. ((2.0)
We now compute the gradient of Jns (0; with respect 0. In the standard policy gradient derivation, the gradient can be swapped with the integral. In this case, the domain of integration depends on thus requiring the Leibniz rule:
V ''
VoJiiss(0;E)
p(z 0)2z
(2.0)
1
<= -:000110000000000000000000111 (4)
V Vep(z 0)
(2.0(
We similarly take the gradient of Equation 3
0 = V
plz(0)dz
=DDo
p(=10)d=
(2,00
(5)
Q (2:0)
Plugging Equation 5 into Equation 4 and rearranging yields
VoJask (0:8)
V Vep(=10) (2-Q1=((2;;)))
(11)
Using the "log-derivative trick" (multiplying by pl Z|) (20))/((()) and using the derivative of loga- rithm) and the definition of conditional expectation yields the final result:
VoJisu
Q (2.0) (i- (2: p(2 log
=E =Ei>p(>10) - (2:
Z Q1=:
The extension the case in which 7 is replaced by scalar reward function R(T). where is generated from parameterized distribution p(T(0), follows the proof in Tamar et al. (2014) without additional adaptation. Thus,
Jrski ( (f; E) (T [(R(T) Rg (0) V log p(T(0) R(T) 2 (,)
where Re is the (1 8) quaattie of the distribution of rewards under the current policy,
APPENDIX C DESCRIPTIONS OF BASELINE ALGORITHMS
Priority queue training. The PQT baseline was implemented using our DSR framework (including situ constraints, hierarchical RNN inputs, and constant optimization), but replacing the risk-secking policy gradient objective with the priority queue training objective proposed in Abolatia et al. (2018 8). PQT works by maintaining a priority queue, defined by the top highest reward samples ever encountered during training, which is updated each training step. They then define their training objective as the average log-likelihood of the samples the priority queue:
10or(0:h) =
1
log p(
k
14

---

Published as a conference paper at ICLR 2021
where & is a hyperparameter controlling the size of the priority queue, and r is the 1thh member of the
priority queue. i.c. the 1t best expression encountered during training. Unlike Jaa(() and Jusk (0:=));)), E), the terms in Jppr(0: (0; k) are not scaled by rewards.
As discussed in [Abolatia et al. 2018), the PQT objective can be combined with a policy gradient objective; however, the authors did not find this to improve results relative to PQT only. For this reason, and to avoid additional hyperparameters, our PQT baseline does not include a policy gradient term.
Vanilla policy gradient. The VPG baseline was implemented using our DSR framework (including in situ constraints, hierarchical RNN inputs, and constant optimization), but replacing the risk-seeking policy gradient objective, Jásk E). with the standard policy gradient objective, Jsa (0). .) In addition, as is common for standard policy gradients, we subtract baseline from the reward function. As in Zoph & Le (017). our baseline is an exponentially- weighted moving average (EWMA) of rewards. This introduces a VPG-specific hyperparameter B3, controlling the degree of weighting decrease. Note that VPG is essentially an ablation of the risk-seeking policy gradient in DSR. However, as it is the critical ablation, we independently tune hyperparameters for the VPG baseline, in contrast to the "No
risk" ablation in Figure 3
Genetic programming. GP-based symbolic regression was implemented using the open-source software package "deap" (Fortin et al. 2012). The initial population of expressions is generated using the "full" method (Koza) 1992) with depth randomly selected between dmin and dmax. The selection operator is defined by deterministic tournament selection, in which the expression with the best fitness among k randomly selected expressions is chosen. The crossover operator is defined by swapping random subtrees between two expressions. The point mutation operator is defined by replacing a random subtree with a new subtree initialized using the "full" method with depth randomly selected between dmin and dmax-
Constraints were implemented in GP in post hoc manner. Each time a new expression is generated via mutation or crossover, all constraints are checked. If any constraints are violated, the new expression is rejected and the evolutionary operation is simply reverted (Fortin et al., 2012 Craenen al. 200T). Like the in situ constraints used in DSR, this process reduces the search space; however, unlike in situ constraints, this rejection-basee process introduces a tradeoff between search space and sample diversity. Thus, we performed all GP hyperparameter grid search both with and without post hoc constraints (using the same constraints as DSR and the RNN-based baselines). Using the best hyperparameters found for each mode, we repeated the full recovery rate experiments on the Nguyen benehmarks, and ultimately selected the best performing mode. We found that constraints actually hindered performance: GP recovery rate with constraints fell from 60.1% % (in Table [ to 32.3% across the Nguyen benchmarks. Thus, GP results herein do not include constraints other than constraining the maximum length to 30 (to prevent excessively long expressions with a low chance of recovery) and the maximum number of constants to 3 (otherwise, runs using constant optimization can become prohibitively expensive) as in DSR.
Eureqa. Considered be the gold standard for symbolic regression (Udreseu & Tegmark 2020), Eureqa is based on the Ga-based approach proposed by /Schmidt& & Lipson (2009). Eureqa experiments were run using the DataRobot platform (www.datarobol
Eureqa allows the user select which tokens to include, and thus experiments used the same library as DSR, as specified by Table 2 in Appendix(] [ While the number of expression evaluations is not tunable or made transparent to the user, Eureqa allows the user to select among several runtime modes. For all experiments, we used the longest-runnin mode, Long Running Search (10, 000 Generat inns) ). After training, Eureqa outputs the best found expression based on mean-square error (MSE), as well as the Pareto front (using a proprietary complexity measure). For noise experiments, to mimic the other baselines, we checked for exact symbolic equivalence for all solutions along the Pareto front.
Wolfram. Wolfram results were obtained using the Mathematica built-in function FindFo rmuuuam la which combines Markov chain Monte Carlo sampling methods with nonlinear regression (Fortuna 2013). The function allows the user to specify the set of tokens to include. However, subtraction and division operators are not allowed. The use of real-valued constants is an integral part of the method and cannot be deactivated. Therefore, the Wolfram experiments used a slightly different library than the other algorithms. Specifically, we used {+, pow, sin, cos, exp, log, 2). where pow(a, b) = a².
15

---

Published as a conference paper at ICLR 2021
Note that b can be obtained via + (-1.0) x and a b can be obtained via pow(b, - -1.0);; thus, all benchmarks can still be recovered using this modified library.
Notable exclusions. The baselines selected for this work are suitable for searching the space of tractable mathematical expressions to discover exact symbolic expressions, and can accommodate numerical constants within expressions. There several well-known algorithms that do not exhibit these features, and thus we do not include as baselines:
Geometric-semante genetic programming (GSGP): Originally proposed by Moraglio et al. (2012). GSGP is a GP-based approach that directly searches the space of the semantics of program. They do this by constructing geometric semantic operators, which promote a strict increase in fitness per generation. However, this desirable feature comes with the cost of the size of the expressions growing exponentially. The resulting expressions routinely on the order trillions of nodes (Pawlak) 2016). This essentially forfeits the ability to exactly recover expressions, and severely limits the usefulness of the technique for interpretabilityy the primary motivation for symbolic regression.
Equation learner (EQL*): As discussed in Related Work, Sahoo al. 2018) perform symbolic regression via neural networks whose activation functions are elementary operators. However, the approach precludes the ability to recovery many simple classes of expressions, such as those involving exponent, logarithm, roots, or division within unary operators. Further, the resulting expressions can be very long, including real-valued shifts and scaled for each operator.
AI Feynman (AIF): As discussed in Related Work, Udrescu & Tegmark (2020) develop tool to recursively simplify symbolic regression problems into smaller sub-problems; solving each sub-problem then requires an inner search algorithm. The authors demonstrate that many physies equations can be reduced to sub-problems that are solvable using simple inner search algorithm, e.g. polynomial fit or small brute-force search. However, more challenging sub-problems- procuuii-incluuding many of the benchmark expressions in this work (namely, those with numerical constants, or those which are pooo--sparable---ae still require an underlying symbolic regression method to conduct the discrete search. Thus, AIF can be used during pre-processing as a siaieemm-imlliiicattii symbolic regression algorithm. Here, DSR and the baselines tested focus on the underlying discrete search problem, c.g. after any problem simplification steps have been applied.
step, then combined with any
APPENDIX D ADDITIONAL EXPERIMENT DETAILS
Detailed description of the symbolic regression benchmarks. Details of the benchmark symbolic regression problems are shown in Table 2 Note that benchmarks without the "const" token in the library do note use a constant optimizer (line in Algorithm [ detailed in 3 in Appendix (A). Benchmarks without real-valued constants Appendix(E) [E), thus recovery is defined by exact symbolic equivalence. Note that the square root and power operators with the given function set, e.g. Nguyen-8 can still be recovered via exp( log(z)) and Nguyen-11 can still be recovered via exp(y log(z)).
be recovered exactly (except Neat-6, discussed in
part of the function set; however, such benchmarks can still be recovered
While all algorithms always produce syntactically valid expressions, they do not always produce semantically meaningful expressions. For example, the expression log(z) is not semantically mean- ingful for non con-positive values of z. Thus, for all commercial algorithms, any expression that produces an overflow or other floating-poin error on the input domain is deemed "invalid" and receives a reward of 0, which corresponds to infinite error. Fureqa and Wolfram follow their own strategy to deal with this issue.
Hyperparameter selection. Hyperparameters were tuned by performing grid search on benchmarks 7 and Nguyen-10. For each hyperparameter combination, we performed 10 independent training runs of the algorithm for 1M total expression evaluations. We selected the hyperparameter combination with the highest average recovery rate, with ties broken by lowest average NRMSE. For all algorithms, the best found hyperparameters were used for all experiments and all benchmark expressions.
16

---

Published as a conference paper at ICLR 2021
Table 2: Benchmark symbolic regression problem specifications. Input variables are denoted by I and/or U(a, denotes e random points uniformly sampled between and for each input variable; training and test datasets use different random seeds. E(a, b, c) denotes c points evenly spaced between and b for each input variable; training and test datasets use the same points (except Neat-6, which uses B(1, 120, 120) as test data, and the Jin tests, which use test data). To simplify notation, libraries are defined relative to a "base" library
U U(-3.3,30) Lo = (+, x, +, sin, cos, exp, log, z}. Placcholder operands are denoted by o, c.g. corresponds to the square operator.
Expression x²
Dataset
Library
Nguyen-1
U,-1, 20)
La
Nguyen-2
1,-1,1,20)
Ca Cg Cg La Co
Nguyen-3
+ x. w² + x² +
U(-1,1,
Nguyen-4
+ x + +
1,-1,1,20)
Nguyen-5
sin(a?) noa(2))
1(-1,1,20)
Nguyen-6
sin(m) san
U(-1,
Nguyen-7
log(z + log( +)
U(0, 2,
Ca
Nguyen-8
V sin(x)) sin(y")
D(0,4,20
Cg
Nguyen-9
10)
(v)
Nguyen-10
2 2sin(x) (o)(
1,0,
Co (y)
Nguyen-11
U(0, 2,
Cg (v)
Nguyen
2,0)
[v]
Nguyen-2'
Au 3x" 2x
U(-1,1,20)
Lo Ca
Nguyen-5'
cos(2) sin( (
U(-
Nguyen 8
V
U(0, 4, 20)
Co Ca Zo [const)
Nguyen-8"
V(0, 4,0,4,20)
Nguyen-1
3.39x 2.12a
1.78x
U(-1,1,20)
Nguyen:S
sin(2 cus(ar) 0.75
U(-1,1,20)
Co [const]
Nguyen 7
log(z 1.4) log(a? (a: 1.3)
D(0, 20)
Co (const)
Nguyen-8
v1.23z
U(0, 20)
Co (const)
Nguyen 10
Sin(115))ccus00.33 (0.5g)
U(0, 2,20)
Co {v. const
Grammar VAE-1
+ sin(z*)
E(-10, 10, 10")
Lo (- 000, log) 2.3
1.3z + 0.5y" 1.7y
U(-3,3, 100]
Lo (log)
const
Jin-2
8.0a² + 8.0y 15.0
U(-3.3,
Lo {log)
const
Jin-3
0.2 + 0.5y 1.2g 0.5
U(-3.3,100 100)
La {log)
const
Jin-4
1.5 exp()( + c.0
U[-3,3,100)
Lo {log] U
const
Jin-5
6.0sin(a) cos(g)
U[-3,3, 100)
Co [log]
const)
6
1.35zy + 5.5 1.0)(y 1.0))
U[-3,3,
Co (log]
V. const
Near-1
+
0(-1,1,20)
(1)
Neat-2
6 xt 3 x
1,-1,
Co (1)
Neat-3
sin(=²) coo(2) -
U(-1,1,20)
Co U (1)
Neat-4
log(z + 1) log( +1)
U(0.2,20)
Co (1)
5
2sin(a) cos(g)
U,-1,1,
(v)
Near-6
El 2 2.1 cos (9.8z) sin(1
E(1, 50,
{+, x,
Neat-7
B(-50, 50, 10")
Le U (tan, tauh,
(e 1)a
Neat-8
U[0.3, 4, 100)
(+,- X, 4, exp.e u (y)
1,24x-2.55 142 14x-4
Neat-9
2(25,,,21)
Conveniently, the three RNN-based algorithms (DSR, PQT, and VPG) each have unique hy- perparameter: risk factor = for DSR, priority queue size k for PQT, and EWMA coefficient | for VPG. The remaining three hyperparameters considered (batch size, learning rate, and entropy weight) are present for each algorithm. For these three algorithms, the space of hyperparameters considered was batch size E (250, 500, 1000), learning rate {0.0003, 0. 0.0005, 0.001), }. and entropy weight Au E {0.01, 0.05, 0.1}. For the algorithm-specinn
hyperparameters, we considered risk factor s E (0.05, 00 0.15) for DSR, priority queue size & (5, 10, 20} for PQT, and EWMA coefficient / {0.1, 0.25, 0.5) for VPG. Thus, these algorithms each considered 81 hyperparameter combinations. The final tuned hyperparameters are listed in Table 3 (Note that hyperparameters tuned independently for each algorithm; identical values across algorithms are incidental.)
For GP, the space of hyperparameters considered was population size E {100, 250, 500, 1000), tourna- ment size {2,3. 3. 10). mutation probability E {0.01, 0.03, 0.05, 0.10, 0.15). crossover probability (0.25, 0.50, 0.75, 0.90, 0.95). and post hoc constraints E { (TRUE, FALSE) (800 combinations).
The final tuned hyperparameters are listed in Table
As commercial software, Eureqa and Wolfram do not expose tunable hyperparameters. Wolfram al- lows the user to select among three "performance goal" modes (Error, Score, and Compl ex exity)
17

---

Published as a conference paper at ICLR 2021
and two "specificity goal" modes (Low and High). We tested all 6 combinations for 100 random seeds. Wolfram is not applicable to Nguyen-10, and we found that recovery rate for Nguyen-7 was 0% for all combinations; thus, the best combination of modalities was selected based on benchmarks Nguyen-5 and Nguyen-6 instead. The best results were obtained using Error for performance and H gh for specificity.
Table 3: Tuned hyperparameters for RNN-based algorithms.
Parameter name
Symbol
DSR
PQT
VPG
Batch size
N
1000
1000
1000
Learning rate
a
0.0005
0.0005
0.0001
Entropy coefficient
0.005
0.005
0.005
Risk factor
E
0.05
Priority queue size
k
10
EWMA coefficient
0.25
Table 4: Tuned hyperparameters for GP.
Parameter
Value
Population size
1,000
Fitness function
NRMSE
Initialization method
Full
Selection type
Tournament
Tournament size (Ar)
Mutation probability
0.05
Crossover probability
0.95
Post hoc constraints
FALSE
Minimum subtree depth (dain)
0
Maximum subtree depth (daax)
Additional details for constant optimization. The same constant optimizer used for all non- commercial software algorithms: DSR, PQT, VPG, and GP. To optimize constants, the values for constant placeholder tokens are optimized against the reward (or fitness) function using BFGS with an initial guess of 1.0 for each constant. Before training, we ensured that all benchmarks with constants do not get stuck in a poor local optimum when optimizing with BFGS and the candidate functional form is correct. Since numerical constants can only be recovered up to floating-poin precision, for benchmarks with constants we determined recovery by manually inspecting the functional form for symbolic correctness. Since constant optimization is a computational bottleneck, we limited each expression to three constants for all experiments with the "const" token, and ran experiments for 1M expression evaluations instead of 2M. Eureqa and Wolfram follow their own strategy to learn real-valued constants.
Additional details for ablation studies. In Figure "No risk" denotes using the standard policy gradient (with baseline) instead of the risk-sceking policy gradient, equivalent to E = 1. Since the standard policy gradient typically includes baseline term, used exponentially-weighhe moving average (weight 0.5) of the average reward of the batch. The difference between this ablation and the VPG baseline is that VPG hyperparameters were tuned independently, whereas this ablation uses the DSR tuned hyperparameters. "No entropy" denotes no entropy bonus, equivalent to Xu O. "No constraints" denotes no constraints precluding nested trigonometric operators, inverse unary operators, minimum length, or maximum length. Instead, if the maximum length of 30 tokens reached, the expression is appended with x until complete. "No risk/ent." denotes combining ablations for "No risk" and "No entropy." "No risk/con." denotes combining ablations for "No risk" and "No constraints." " "No ent./con." denotes combining ablations for "No entropy" and "No constraints." "All ablations" denotes combining ablations for "No risk." "No entropy," and "No constraints."
Pareto front computation. Overfitting in symbolic regression occurs when expression has added complexity to conform to noisy data. For example, high-degree polynomials can fit small datasets with low error, though generalization may be poor. While overfitting is not issue for noiseless experiments (since the global optimal expression has zero error), it is indeed possible for expressions to overfit in the noisy data experiments, especially with small dataset sizes. Thus, for noisy data
18

---

Published as a conference paper at ICLR 2021
experiments, defined recovery exact symbolic equivalence any expression along the Pareto front (defined over the reward-come complexxty plane) at the end of training. To compute the Pareto front, we define simple complexity measure C:
T
C(T) Σ c(T)),
where is the complexity of particular token. In accordance with Eureqa's default parameters, we used token complexities of for +, ; x, input variables, and constants; for +; for sin and cos; and 4 for exp and log.
Note that C' is only used to determine the Pareto front at the end of training; does not affect training. This method of determining recovery based on the Pareto front was used for all algorithms, except that Eureqa uses a proprictary complexity measure.
Computing infrastructure. Experiments were executed on an Intel Xeon E5-2695 v4 equipped with NVIDIA Tesla P100 GPUs, with 32 cores per node, GPUs per node, and 256 GB RAM per node.
APPENDIX E COMPARISONS TO LITERATURE-REPORIT RESULTS
Unfortunately, the majority symbolic regression methods are inconsistent in their experimental setup, i.e. choice of library, benchmarks, and measure of performance. Often, benchmark expressions are hand-selected from several benchmark suites [Tujillo al., 2016) or even hand-crafted for particular study (Kusner et al.) 2017 Sahoo et al., 2018). The criteria for determining "recovery" (also called "hit" or "success") is often defined by an arbitrary error threshold instead of exact symbolic equivalence (White et al. 2013). Further, different works use different performance metrics (e.g. RMSE, mean absolute error) that cannot be converted without access to raw data (Jin et al. 2019).
results. These issues make it difficult to compare methods via direct comparison of literature- To enable such comparisons, we mimic the experimental setup and measures of performance for several works. The alternative approach would be to evaluate these methods using experimental setup and measures of performance; however, not all methods have available source code or suffi- cient implementation details for reproducibility. Below, recapitulate the experiments of several additional works, and directly compare literature-reporteed results to our experiments using DSR. We note that even these comparisons are imperfect for for benchmarks with randomly sampled points, as the random number generator and seed used to generate the training and test data are typically not disclosed.
Semantics- based symbolic regression. Huynh et al. (2016) is an exceptional case in that they also report recovery as exact symbolic recovery, using an overlapping set of benchmarks: Nguyen-1 through Nguyen-10. In this work, the authors present Semantic- based Symbolic Regression (SSR), Ga-based approach that considers the semantics of each expression tree; for example, they remove duplicate trees with different tree structures but identical semantics, e.g. ((((x- 1) + (x/zx) -) + equivalent to y. In Table we compare DSR's recovery to SSR's reported values. To replicate their experimental setup, report DSR recovery rate for 30 independent training runs. DSR's ' average recovery rate is more than twice as high as SSR.
Table 5: Recovery rate comparison of DSR and literature-reported values from SSR (Huynh et al. 2016 A bold value represents statistical significance (p < 10 10-3) across all benchmarks.
Benchmark
Expression
DSR
Nguyen-1
+ x² +x
100%
100%
Nguyen-2
4
+ 3 + x² + r x. + I² + t x² + +x
100%
10%
Nguyen-3
+
100%
70%
N
6 " +
100%
30%
5
sin(x² cos(x) -
77%
10%
Nguyen-6
sin(x) sin(x 22)
100%
73%
Nguyen-7
log(æ + 1) + log(2² +1)
33%
3%
Nguyen-8
VI sin(x) + sin(y?) 2sin(2r) cos( (n)
97%
3%
Nguyen-9
100%
57%
Nguyen-10
100%
93%
Average
90.7%
45.0%
19

---

Published as a conference paper at ICLR 2021
GrammarVAl VAE. As mentioned in Related Work, Kusner et al. (201)) develop generative model for discrete objects that adhere to pre-specified grammar, then optimize them latent space. To mimic their experimental setup. run 10 independent trials using the performance metric log(1 + MSE) and the Grammar VAF benchmark described in Table 2 Using the best solution found in each trial, DSR's average performance metric is 0.0105 + 0.0149, whereas GrammarVAR reported average is orders of magnitude larger: 3.47 + 0.24. In Table 6 as in Kusner et al. 2017). we compare the top expressions reported in GrammarVAL In addition to exactly recovering the expression, DSR'T top 3 expressions are all superior to those found by GrammarVAE Further, when repeating for 100 trials, we find that DSR exactly recovers the benchmark expression in 19% of runs, whereas GrammarVAB never exactly recovers the ground truth expression.
the top expressions found using DSR across all 10 trials.
Table 6: Comparison of DSR and literaturereeported values from Grammar VAE Kusner et al. 2017
Algorithm
Rank
Expression
log(1 + MSE)
GrammarVAN
1
sin(3) + sin(x")
0,04
GrammarVAs VAE
2
1 + sin xx
0.10
GrammarVAIF
3
1 + x + sin(x") (x²
0.37
DSR
1
+ sin( (x*)
sin(1.5)
DSR
2
sin(x
7.0 10
DSR
3
sin (1) + x + (in
3.8 x 10
Bayesian Symbolic Regression. In Jin et al. (2019), the authors propose Bayesian symbolic regres- sion (BSR), a Bayesian framework for symbolic regression. In BSR, several carefully designed prior distributions are used to encode domain knowledge like preference basis functions or tree structure. The resulting posterior distributions efficiently sampled using Markov chain Monte Carlo techniques. To mimic their experiments, we run 50 independent trials using the performance metric RMSE and the Jin benchmarks deseribed in Table Table 7 shows the average RMSE on the test data for DSR and BSR across the benchmarks used in Jin et al. (2019). DSR outperforms BSR on all six benchmarks.
Table 7: RMSE comparison of DSR and litcrature-reported values from BSR (Jin et al.
2019
Benchmark
Expression
DSR
BSR
Jin-1
2.5x 1.3x + 0.5y² 1.7y
0.46 + 0.41
2.04 +: 3.27
Jin-2
8.0z² + 8.0y² 15.0
0
6.84 10.10
Jin-3
0.2* + 0.5y² 1.2y 0.5x
5.2 x 10
4 3.5 X 10
Jin-4
1.5 exp(x) + 5.0 cos
1.36
10
- 8.91 X 10
0.16 0.62
Jin-5
6.0 (.0sin(z)) cos(g)
0
0.66 + 1.13
Jin-6
1.35ry + (i5 (y0)
1.0)) Average
2.23 + 0.94
4.63 + 0.62
0.45 + 0.23
2.42 + 2.66
Neat-GP. In Trujillo et al. (2016), the authors propose Neat-GP, a genetic programming approach that uses the NEAT (NeuroEvolution of Augmenting Topologies) algorithm to address the problem of uncontrolled program size growth standard GP approaches. The authors report median values of RMSE on the test data for 30 independent trials using the Neat benchmark set described in Table 2
We replicate their experimental setup using DSR and compare to their reported values Table DSR outperforms Neat-GP seven of the nine benchmarks.
The task in 6 is to find closed-form approximations the harmonic series H,, = =E= 1. As additional experiment, expressions which are remarkably accurate approximations of H,; when extrapolating to all n e N, the error of the best discovered expression is less than 0.000001% In particular, DSR discovered the expression:
repeated DSR on Neat-6 with (log, const) added to L. DSR discovered
Hn 7 log(n) + 2n 11.3776 12.
+ 0.327981
where y Rs 0.57721 is the Mascherrni constant. This is variation of the formula H,, s + log(n) + 1/(2n) - 1/12n found by Leonhard Euler 1755 (Bromwich] 1908). Remarkably, the constant naturally emerged in our recovered expression.
20

---

Published as a conference paper at ICLR 2021
Table 8: Comparison of median values of RMSE for DSR and literature-reported values from Neat-GP
(Trjiillo et al. 2016).
Benchmark
Expression
DSR
Neat-GP
Neat-1
x x²
+x
0
0.0779
Neat-2
2 + x
0.0576
Neat-3
sin(x²) cos(2))
0.0041
0.0065
Neat-4
log(z
12g(2" +1)
0.0189
0.0253
Neat-5
2 (z) cos (os
0.0023
Neat-6
Σ k=
0.2378
0.2855
Neat-7
2.1 (9.8z) (1.3y)
1.0606
1.0541
Neat-8
0.1076
0,1498
122-2-2.662 1+z
Neat-9
+ 1+y
0.1511
0.1202
Average
0.1756
0.1977
APPENDIX F ADDITIONAL RESULTS ON NGUYEN BENCHMARKS
Variations of Nguyen benchmarks. We introduce two sets of variations the Nguyen benchmarks. In the first set, alter existing benchmarks to make them more challenging. In the second set, introduce real-valued constants and add the "const" token to the library to demonstrate the constant optimizer. We report recovery rate for DSR and all five baselines in Table DSR significantly outperforms all baselines in the first set, and achieves 100% recover in the second set.
Table Recovery rate comparison DSR and five baselines on variations of the Nguyen symbolic regression benchmark suite, A bold value represents statistical significance (p < 10 10-$) across each set of benchmarks.
Benchmark
Expression
DSR
PQT
VPG
GP
Eureqa
Wolfram
Nguyen-2'
4x + 3x 2x
1 +
96%
92%
0%
52%
0%
100%
N
sin(x²)
cos(x:) 2
87%
25%
0%
4%
1%
0%
Nguyen-8'
50%
44%
1%
2%
0%
0%
Nguyen-8"
3%
1%
0%
1%
1%
0%
Average
59.0%
40.5%
0.2%
14.8%
0.5%
25.0%
Nguyen-1
3.39x + 2. 2.12x² + 1.78r
100%
100%
100%
100%
100%
100%
N
sin(x² cos(x) 0.75
100%
100%
100%
70%
100%
1%
Nguyen-7
log(x + 1.4) + 1xg(x² 1.3)
100%
90%
100%
0%
0%
0%
N
VI.23r
100%
100%
100%
20%
30%
Nguyen-10
sin(1.5xr cos (0.5g)
100%
100%
0%
10%
Average
100.0%
98.0%
88.0%
38.0%
48.0%
NRMSE comparisons. Tables and show recovery rates for all algorithms across all benchmarks. In Table I0 we show the analogous results for NRMSE on the test data. DSR outperforms all algorithms across each benchmark set, except Eureqa on the original Nguyen benchmark set. Here, Eureqa's 's lower average NRMSE is attributed to its low error for Nguyen-12. However, this is a result of not being able constrain or limit the complexity of expressions produced by Eureqa. In particular, Eureqa routinely identified expressions of length ~00 for Nguyen-12 12, whereas all non-commercial baselines were limited to length 30 ensure that expressions are tractable.
We also observe that for the few expressions with low or zero recovery rate (e.g. Nguyen-7 and N2uyen-12), GP sometimes exhibits lower NRMSE. One explanation is that GP is more prone to overfitting the expression to the dataset. As evolutionary approach, GP directly modifies the previous generation's expressions, allowing it to make small "corrections" that decrease error each generation even if the functional form is far from correct. In contrast, in DSR the RNN "rewrites" each expression from seratch each iteration after learning from gradient update, making it less prone to overfitting.
Training curves. Figures 6 and 7 show the reward (11++NRMME) and recovery rate, respectively, function of total expressions evaluated during training. Commercial software algorithms (Fureqa
21

---

Published as a conference paper at ICLR 2021
Table 10: Comparison of NRMSE on the test data for DSR and five baselines on original and variations of the Nguyen symbolic regression benchmark suite.
Renchmark
DSR
PQT
VPG
GP
Farega
Wolfram
Nguyen-1 1
0.000 1 0.000
0.000
1: 0.000
0.001 1 0.004
0.000 1 0.000
0.000 1: 0.000
0.000 1: 0.000
Nguyen-2
0.000 +: 0.000
0.001
0.005
0.021
0.051
0.001 0.003
0.000 1: 0.000
0.000 + 0.000
Nguyen-3
0.000 = 0.000
0.004 = 0.010
0.051 + 0.045
0.000 = 0.000
0.001 = 0.005
0.000 + 0.000
N
0.000 1: 0.000
0.001 1: 0.004
0.041 0.02
0.000 1 0.000
0.003 1: 0.006
0.000 1: 0.000
Nguyen-5
0.015 = 0.026
0.007 -007=00 0.014
0.022
0.005 0.011
0.001 - 0.003
0.502 + 0.398
Nguyen-6
0.000 = 0.000
0.000 0.003
0.000
0.001 = 0.002
0.000 - 0.000
0.609 + 0.059
Nguyen-7
0.005
0.007
0.001
0.003 0.001
0
Nguyen 8
0.003 1: 0.014
0.048 1. 0.037
0.125
0.072 1 0.081
0.075 0.115
0.006 1 0.034
Nguyen-9
0.000 + 0.000
0.000 0.000
0.000
0.000 0.000
0.000
Nguyen-10
0.000 = 0.000
0.004
0.001 0.001
0.006 = 0.018
0.006
Nguyen- II
0.000 + 0.000
0.000 + 0.000
0.081 = 0.15
0.000
Nguyen-12
0.202 0.060
0.212
0.249 +
0.139 0.041
0.068 #
Average
0.019 = 0.059
0.023 + 0.062
0.018 + 0.082
0.026 = 0.065
0.013 + 0.016
Nguyea-2'
0.001 = 0.006
0.003 + 0.010
0,190 + 0.071
0,023 = 0.101
0.018 + 0.020
0.000 + 0.000
Nguyen-5'
0,002 + 0.004
0.022 0.033
0.221 0.160
0.022 + 0.017
0.011 + 0.012
0.901 + 0.297
Nguyca-8'
0.016 0.121
0.083 086
0.322 0.269
0.173 0.230
0.113 - 0.123
0.218 + 0.106
Nguyen.s"
0.054 + 0.028
0.041 + 0.027
0.089 + 0.023
0.047 + 0.068
0.029 + 0.038
0.041 + 0.089
Average
0.026 = 0.066
0.036 0.097
0.201 0.179
0.000 = 0.144
0.042 = 0.077
0,290 - 0.443
Nguyen-1"
0.000 i: 0.000
0.000 1 0.000
0.000 0.000
0.000 1 0.000
0.000 1 0.000
0.000 1 0.000
N5uyen-5
0.000 + 0.000
0.000 0.000
0.000 0.000
0.000 0.000
0.000 4: 0.000
0.714 + 0.449
Nguyen-7e
0.000 = 0.000
0.001 + 0.002
0.000 0.000
0.011 = 0.031
0.004 = 0.016
0.043 + 0.000
Nguyen-8*
0.000 1: 0.000
0.000 0.000
0.000 0.000
0.005 0.004
0.002 1: 0.003
0.045 1 0.073
Nguyen-10"
0.000 + 0.000
0.000 + 0.000
0.004 0.006
0.002 - 0.001
0.000 + 0.000
Average
0.000 = 0.000
0.000 0.001
0.001 = 0.003
0.004 = 0.014
0.001 = 0.007
and Wolfram) only provide the final solution, so training curves not possible. Note that the recovery rate values in Table I correspond to the final point on each curve in Figure
Nguyen-1
Nguyen-2
Nguyen-3
Nguyen-4
1.0
Reward
0.9
DSR
PQT
0.8
VPG
Nguyen-5
Nguyen-6
Nguyen-7
Nguyen-8
1.0
Reward
0.9
Nguyen-9
Nguyen-10
Nguyen-11
Nguyen-12
1.0
Reward
0.9
0.8
500
1000
1500
2000
500 1000
1500
2000
500 1000
1500
2000
500
1000
1500
2000
Training step
Training step
Training step
Training step
Figure 6: Reward training curves for DSR, PQT, VPG, and GP for the Nguyen benchmarks. Each curve shows the best reward (1/(1 + NRMSE)) found so far as function of expressions evaluated, averaged across all training runs. A value of 1.0 denotes that all training runs recovered the correct expression. Error bands represent standard deviation.
Distributions of rewards during training. Figure 2 compares the performance of the risk-secking policy gradient and standard policy gradient for Nguyen-8 by showing the distributions of rewards during training. Analogous plots for all 12 Nguyen benchmarks are shown in Figures and
22

---

Published as a conference paper at ICLR 2021
Nguyen-1
Nguyen-2
Nguyen-3
Nguyen-4
1.00
Recovery
0.75
DSR
PQT
0.25
VPG
Nguyen-5
Nguyen-6
Nguyen-7
Nguyen-8
1.00
Recovery
0.75
0.50
0.25
Nguyen-9
Nguyen-10
Nguyen-11
Nguyen-12
1.00
Recovery
0.50
0.25
0.00
500
Training step
Training
Training step
Training
Figure 7: Recovery rate training curves for DSR, PQT, VPG, and GP for the Nguyen benchmarks. Each curve shows the fraction of independent training runs that correctly recovered the benchmark expression as a function expressions evaluated. value of 100% denotes that all training runs recovered the correct expression.
APPENDIX G RUNTIMES OF NON-COMMERCIAL ALGORITHMM
All experiments used an identical number total evaluation expressions for each method. Another consideration is the compute time of different methods. In Table 11 report the average runtimes (on a single core) for non-commercial algorithms across the Nguyen benchmark set. We report two runtime metrics: 1) runtime process all 2M expressions, without early stopping if the benchmark expression is recovered: 2) runtime until recovering the benchmark expression. Without early stopping, GP is unsurprisingly the fastest method, it does not require neural network training. However, when considering early stopping, DSR is the fastest, largely because it has the highest recovery rate, thus triggering early stopping more often.
Table 11: Comparison of runtimes across the Nguyen benchmark set.
Early stopping
Method
No
Yes
DSR
483.0 = 641.9
PQT
2320.3 + 2
959.6 = 925.9
VPG
2452.3 107.4s
1477.1 1058.7
GP
1375.4 + #34.0
915.5 + 4
23

---

Published as a conference paper at ICLR 2021
Full basc
20-
Top (anch
Nguyen-2
Full batch
Top (tatch
20
12
Standard boir
Ful torm
Tep Barth
caanh
0.
Rewort
Revove
Reppnt
Kewere
10-
Nguyen-3
Miik aiit-nssing nG.
Risk swaking PG.
10-
Niuuye-4
Rakissedin PC.
Ral suoidng PG.
10-
Standu
Stretir No,
10
Standard
Surtlad
Dundy
20
0.0
Nguyen-5
Full bescn
20
15
Top e iteich
26
Nguuynn-6
Full butch
26
Top chaich
32
26-
Ful basch
Tops barch
Full batch
12
Top batch
20
0.0
6.2
0.4
Roward
Revord
Reward
Roward
Nguyen-7
Hix canking Po.
Aikk-ceed eg Ph.
Niguyen-8
takaseding
nal-sssiing PG
Standay
Staddrtd is.
Stndart Po.
Nguuyen
Tull busch
20
Top balsh
28
Nguuynn110
Fut butch
Tea abatch
I6-
Hurdard
Standerd sharch
Fullwach
Tap Eatth
nd borth
10
Revard
Ngu 11
Nguyen-12
Raksuedrp PC.
Ital auaiiin
betch
Standr
Standard
Sandard Top
5.0
Full batch
Tep batch
642
0.8
1.0
60
42
0.4
3.0
0.6
1.0
0.5
0.2
7A
2.6
0.5
1.0
0.2
0.2
0.4
0.6
0.5
Rovard
Reward
Figure 8: Empirical reward distributions for Nguyen benchmarks, with and without risk-secking policy gradients. Each group of four plots corresponds particular benchmark expression. Each curve is a Gaussian kernel density estimate (bandwidth 0.25) of the rewards for a particular training iteration, using either the full batch of expressions (plots labeled "Full batch") or the top fraction of the batch (plots labeled "Top batch"), averaged over all training runs. Black plots were trained using the risk-seeking policy gradient objective. Blue plots were trained using the standard policy gradient objective. Colorbars indicate training step. Triangle markings denote the empirical mean of the distribution at the final training step.
24

---

Published as a conference paper at ICLR 2021
Nguyen-1
Nguyen-2
Nguyen-3
Reward
Full
0.2
Risk-seeking PG
Top
Standard PG Best
0.0
Nguyen-4
Nguyen-5
Nguyen-6
Reward
0.6
0.4
0.2
0.0
Nguyen-7
N
Nguyen-9
1.0
0.8
Reward
0.6
0.4
0.2
0.0
Nguyen-10
Nguyen-11
Nguyen-12
1.0
0.8
Reward
0.6
0.4
0.2
0.0
Expressions
Expressions
Expressions
Figure Training curves for Nguyen benchmarks, with and without risk risk-seeking policy gradients. Dotted curves denote mean reward of the full batch. Dashed curves denote mean reward of the top fraction of the batch. Solid curves denote the best expression found so far. Each curve is averaged over all training runs.
25

---

Published as a conference paper at ICLR 2021
REFERENCES
Daniel A Abolatia, Mohammad Norouzi, Jonathan Shen, Rui Zhao, and Quoc V Le. Neural program synthesis with priority queue training, arXiv preprint arXiv: 00356, 2018.
Thomas John I' Ansoo Bromwich, An introduction to the theory of infinite series. Macmillan and Company, limited, 1908.
BGW Craenen, AE Eiben, and E Marchiori. How to handle constraints with evolutionary algorithms. Practical Handbook Of Genetic Algorithms: Applications, pp. 341-361, 2001.
Félix-Antoine Fortin, François-Michel De Rainville, Marc-André Gardner, Marc Parizeau, and Christian Gagné. Deap: Evolutionary algorithms made easy. Journal of Machine Learning Research, 1301))22111-215, 2012.
Giorgia Fortuna. Automatic formula discovery in the wolfram language. In Wolfram Technol- ogy Conference 2015, Wolfram Research, 2015. URL httttt:ss///bbaeeee
folffrm... com/
Srrocentter/Cofeeeeeeee /9329 9329/
Quang Nhat Huynh, Hemant Kumar Singh, and Tapabrata Ray. Improving symbolic regression through a semantics-driven framework. In 2016 IEEE Symposium Series on Computational Intelligence (SSCI), pp. 1-8. IEEE, 2016.
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian symbolic regression. arXiv preprint 1XXi::990000892 2019.
John R Koza. Genetic Programming: On the Progrumming of Computers by Means of Natural Selection, volume MIT press, 1992.
Matt Kusner, Brooks Paige, and José Miguel Hernández- Grammar variational autoencoder. In Proceedings of the 34th International Conference Machine Kearning-Volunne 70, pp. 1945- 1954. JMLR. org, 2017.
Alberto Moraglio, Krzysztof Krawiec, and Colin G Johnson. Geometric semantic genetic program- ming. International Conference on Parallel Problem Solving from Nature, pp. 21-31. Springer, 2012.
Tomasz P Pawlak. Geometric semantic genetic programming is overkill. In European Conference Genetic Programming, pp. 246-260. Springer, 2016.
Subham S Sahoo, Christoph H I Lampert, and Georg Martius. Learning equations for extrapolation and control. arXiv preprint 18Xiv:1806.0740 07259, 2018.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data, science, 324(5923):81-855 2009.
Aviv Tamar, Yonatan Glassner, and Shie Mannor. Policy gradients beyond expectations: Conditional value-at-risk. arXiv preprint 4rXiv:1404.3822, 2014.
Leonardo Trujillo, Luis Muñoz, Edgar Galván-López, and Sara Silva. neat genetic programming: Controlling bloat naturally. Information Sciences, 333:21-43, 2016.
Silviu-Marian Udrescu and Max Tegmark. Ai feynman: physics-inspired method for symbolic regression. Science Advances, 6(16)eeaay6311, 2020.
David R White, James Medermott, Mauro Castelli, Luca Manzoni, Brian W Goldman, Gabriel Kronberger, Wojciech Jaskowski, Una-May O'Reilly, and Sean Luke. Better gp benchmarks: community survey results and proposals. Genetic Programming and Evolvable Machines, 14(1): 2013.
Barret Zoph and Quoc V Le. Neural architecture search with reinforcement learning. International Conference on Learning Representations, 2017.
26


Reviews:
[{"review_id": "fkXSX0BWQn", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper presents a novel approach to the problem of symbolic regression, where the goal is to learn relationships between variables in the form of mathematical expressions. This is clearly a very relevant task towards constructing explainable AI systems.\n\nThe proposed approach in this paper is based on the generation of mathematical expression with recurrent neural networks, by exploiting background knowledge about the form of the expressions to impose constraints on the generated examples.\n\nThe RNN is trained by maximizing a risk-seeking policy gradient that aims to increase best-case performance. The key idea here is to increase the reward of the top-epsilon fraction of samples from the distribution, without taking into account the samples that fall below such threshold.\n\nThe technique is sound and novel, and it provides a significant contribution to this research area. The positioning of the paper with respect to the state-of-the-art in the field is highly accurate.\n\nA very solid experimental evaluation is carried out on several benchmarks, comparing the proposed approach with state-of-the-art systems for the same task, including commercial software such as Eureqa and Wolfram. The analysis includes an ablation study and experiments conducted with different amounts of training data. The proposed methodology is shown to perform better than all the competitors.\n\nOverall, I consider the paper to be a strong contribution for ICLR.\n\n* In the experiments with different levels of noise in data, why was Gaussian noise added to the dependent variable only, and not also to the independent variables?\n\n- Pag. 4, \"is not allowed.While\" -> \"is not allowed. While\"\n- Pag. 5, \"but in practice has high variance\" -> \"but in practice it has high variance\"\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "W9pGYaYB1Q", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary\n\nThe paper uses RNN trained by risk seeking RL objective to predict mathematical expression that generated the target dataset. In the second step placeholders for constants in sampled expression are optimized again by gradient optimizer to maximize the reward function. The decoder RNN uses heuristics that make it easier to generate valid math expressions.\n\n\nStrong points\n\n+ The paper is clearly written and potentially relevant to larger audience.\n+ The proposed method is relatively easy to implement and it leads to good empirical results.\n\n\nRecommendation\n\nI recommend acceptance of this paper since the proposed technique is highly competitive, easy to implement and well presented.\n\n\nQuestions\n\n* How does compute used by GP to get results in Tab 1 compare to policy gradient based methods? My perspective is that all the techniques should find the correct solution in the limit with enough computational resources (under assumption that they can escape local optima). Therefore making sure that all the methods had roughly similar budgets is important.\n\n* In figure 2E mean reward for standard PG is going up even after 2M training steps, would it be possible to run the experiment longer to see where the max would be once mean reward converges? (Similarly in figure 8 in tasks Nguyen 3,5 and 8 mean of standard GP is still going up.) It is great to see that risk seeking PG converges much faster on these tasks, however knowing what is the limit for standard PG would be also interesting.\n\n* In section 3.1 you say that generated expressions had to be between 4 and 30 symbols long, however tasks 8 and 11 (sqrt(x) and x^y) have shorter descriptions. What am I missing?\n\n\nPossible improvements\n\n* Using the set of problems from AI Feynman would make empirical evaluation stronger. However the same can be said for AI Feynman system and Nguyen dataset.\n\n* Beyond scope of this paper: At the moment a new sequential model is learned for each task and it learns about that task only through the reward function. What if the model can \"see\" the dataset first by reading it through n-dimensional CNN (or RNN) that would produce \"dataset embedding\" that will be later used to condition the sampling RNN. In the same way that image captioning model is conditioned on image embedding. Generating infinite training dataset (with random expressions) should be trivial, that is one advantage over image captioning with limited data.\n\nTypos\n\nstandard standard -> standard", "extra_scores": {}, "extra_text": {}}, {"review_id": "6Gazor0EiJ", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "#### Summary\n\nThis paper presents a symbolic regression algorithm which uses policy gradient to learn a distribution over the space of mathematical expression structures. The distribution is represented by an RNN, and the on-policy sampling is also realised by this RNN. To prune the massive sample space, the authors include several constraints according to some prior domain knowledge. Furthermore, instead of optimising the RNN directly with policy gradient, the proposed algorithm also makes use of the risk-constrained method that emphasises the risks above a given percentile criteria. The proposed approach has achieved the best performance on several benchmark symbolic regression datasets.\n\n#### Pros\n\n+ This paper is motivated by an interesting problem. I like the idea of using domain knowledge to introduce hierarchical constraints, which seems to be effective in this task.\n+ The authors have done extensive experiments, the presented algorithm performs well on the benchmark dataset.\n+ This paper has covered a wide range of related work.\n\n#### Cons\n\n- The main contribution of this paper is introducing the heuristic constraints to prune the expression search space and applying existing reinforcement learning techniques (e.g., using RNN to represent the policy of sampling a grammar; the risk-constrained reinforcement learning) to the symbolic regression task. It looks more incremental rather than a novel contribution.\n- The number of experiments in the main article is small, and many experimental details are hidden in the supplementary, I think the authors should re-organise the paper and prompt some of the results in the appendix to the main article.\n- One of the biggest problem for the on-policy reinforcement learning is its time complexity since it requires a large sample for estimating the expectation of risk. Considering that symbolic regression is a hard problem with highly complicated underlying distribution over the expression structures, I wonder how efficient is the proposed approach comparing to the other baseline methods. I think the authors should include the results of running time of all the compared approaches.\n- The experimental setting is questionable for some compared methods. In Appendix E, the authors report the results of several SOTA symbolic regression approaches. However, as the authors stated, the experimental results are different because the choice of datasets and primitive functions are different. I have noticed that this is a re-submission from ICLR 2020, and reviewers were pointing out the problem of lacking comparison one year ago. I think one year should be enough for carrying out experiments of all these SOTA approaches under the same experimental setting. Still, there is no such result in the current submission.\n\n#### Recommendation\n\nI think this paper proposes a nice approach to tackle the symbolic regression problem. Although the presented algorithm seems to be a combination of existing approaches, its performance on benchmark data sets looks quite good. However, the experiments are still weak because this paper lacks a fair comparison (e.g., using same library and the same set of problems), and my biggest concern of the proposed approach is its time complexity. I think the authors should address these problems to improve this work.\n\n--------------------------------------------------------------------------\n#### Update after rebuttal\nThe authors have explained the reasons for the comparing experiments in Appendix E and updated the result of time complexity, which I think worth including in this paper.", "extra_scores": {}, "extra_text": {}}, {"review_id": "-6JSOb3_4rr", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 9, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors introduce a novel method for inferring a simple algebraic expression for an output in terms of an input. The method outputs a sequence tokens of the algebraic expression as represented by the pre-order traversal of an expression tree. For each token of the sequence, a recurrent neural network outputs a probability distribution on possible tokens, from which the token is sampled. Some notable contributions:\n\n1. The authors introduce a novel reinforcement learning objective that optimizes the quality of the best examples instead of optimizing the average quality.\n\n2. At each point in the sequence, their recurrent neural network takes as input the sibling and parent nodes of the expression tree, instead of the previous token in the sequence.\n\n3. Simple constraints can be introduced in their approach by zeroing out the sampling probability for tokens that don't meet the constraints. (For example, they disallow redundant expressions like log(exp(x))).\n\nThe authors do a comprehensive comparison with other methods (in which their work performs very favorably) as well as an ablation analysis to confirm that each of their innovations is helpful (which I was very pleased to see).\n\nWell done.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an approach of generating mathematical expressions with a recurrent neural network, which is trained with risk-seeking policy gradient to maximize the quality of best examples rather than average examples.  The proposed approach also enables easily incorporating domain knowledge or constraints to avoid illegal or redundant expressions.  In extensive experiments, the proposed method is shown to significantly outperform strong baselines, including commercial software.  All of the reviewers find the work interesting and relevant, and there are no major concerns or issues after discussion.  The topic is also of interest to a wide range of audience in the ICLR community.
Decision: Accept (Oral)