Year: 2019

Abstract:
Discovering the underlying mathematical expressions describing a dataset is 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 are lacking. We propose a framework that combines deep learning with symbolic regression via a simple idea: use a large model to search the space of small models. More specifically, we use a recurrent neural network to emit a distribution over tractable mathematical expressions, and employ reinforcement learning to train the network to generate better-fitting expressions. Our algorithm significantly outperforms standard genetic programming-based 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 a priori constraints in situ.

heading: 1 INTRODUCTION
text: Understanding the mathematical relationships among variables in a physical system is an integral component of the scientific process. Symbolic regression aims to identify these relationships by searching over the space of tractable mathematical expressions to best fit a dataset. Specifically, given a dataset of (X, y) pairs, where X ∈ Rn and y ∈ R, symbolic regression aims to identify a function f(X) : Rn → R that minimizes a distance metric D(y, f(X)) between real and predicted values. That is, symbolic regression seeks to find the optimal f? = argminf D (y, f(X)), where the functional form of f is a tractable expression.
The resulting expression f? may be readily interpretable 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. However, 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 an extremely challenging machine learning problem.
Given the large and combinatorial search space, traditional approaches to symbolic regression typically utilize evolutionary algorithms, especially genetic programming (GP) (Koza, 1992; Bäck et al., 2018). In GP-based symbolic regression, a 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., 2015). However, deep learning may seem incongruous with or 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 while entirely bypassing the need to interpret a network.
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 parameters using a policy gradient algorithm. As training proceeds, the RNN adjusts the likelihood of an expression relative to its reward, assigning higher probabilities to better fitting expressions.
We demonstrate that DSR outperforms a standard GP implementation in its ability to recover exact symbolic expressions from data, both with and without added noise. We summarize our contributions as follows: 1) a novel method for solving symbolic regression that outperforms standard GP, 2) an autoregressive generative modeling framework for optimizing hierarchical, variable-length objects, 3) a framework that accommodates in situ constraints, and 4) a novel risk-seeking strategy that optimizes for best-case performance.

heading: 2 RELATED WORK
text: Symbolic regression. Symbolic regression has a long history of evolutionary strategies, especially GP (Koza, 1992; Bäck et al., 2018; Uy et al., 2011). Among non-evolutionary approaches, the recent AI Feynman algorithm (Udrescu & Tegmark, 2019) is a multi-staged approach to symbolic regression leveraging the observation that physical equations often exhibit simplifying properties like multiplicative separability and translational symmetry. The algorithm identifies and exploits such properties to recursively define simplified sub-problems that can eventually be solved using simple techniques like a polynomial fit or small brute force search. Brunton et al. (2016) develop a sparse regression approach to recover nonlinear dynamics equations from data; however, their search space is limited to linear combinations of a library of basis functions.
AutoML. Our framework has many parallels to a body of works within automated machine learning (AutoML) that use an autoregressive RNN to define a distribution over discrete objects and use reinforcement learning to optimize this distribution under a black-box performance metric (Zoph & Le, 2017; Ramachandran et al., 2017; Bello et al., 2017). The key methodological difference to our framework is that these works optimize objects that are both sequential and fixed length. For example, in neural architecture search (Zoph & Le, 2017), an RNN searches the space of neural network architectures, which are 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, a major contribution of our framework is defining a search space that is both inherently hierarchical and variable length.
The most similar AutoML work searches for neural network activation functions (Ramachandran et al., 2017). While the space of activation functions is hierarchical in nature, the authors (rightfully) constrain this space substantially by positing a functional unit that is repeated sequentially, thus restricting their search space back to a fixed-length sequence. This constraint is well-justified for learning activation functions, which tend to exhibit similar hierarchical structures. However, a repeating-unit constraint is not practical for symbolic regression because the ground truth expression may have arbitrary structure.
Autoregressive models. The RNN-based 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;b) in addition to the AutoML 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 (You et al., 2018). In principle, we could have constrained GraphRNN to define the distribution over expressions, since trees are a special case of graphs. However, GraphRNN constructs
graphs breadth-first, whereas expressions are more naturally represented using depth-first traversals (Li et al., 2005). 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 a node’s value determines its number of children (e.g. cosine is a unary node).

heading: 3 METHODS
text: Our overall approach involves representing mathematical expressions by the pre-order traversals of their corresponding symbolic expression trees, developing an autoregressive model to generate expression trees under a pre-specified set of constraints, and using reinforcement learning to train the model to generate better-fitting expressions.

heading: 3.1 GENERATING EXPRESSIONS WITH A RECURRENT NEURAL NETWORK
text: We leverage the fact that algebraic expressions can be represented using symbolic expression trees, a type of binary tree in which nodes map to mathematical operators, input variables, or constants. Operators are internal nodes and may be unary (e.g. sine) or binary (e.g. multiply). Input variables and constants are terminal nodes. We encode an expression τ by the pre-order traversal (i.e. depth-first, then left-to-right) of its corresponding expression tree.1 We denote the ith node in the traversal as τi and the length of the traversal as |τ | = T . Each node has a value within a given library L of possible node values or “tokens,” e.g. {+,−,×,÷, sin, cos, x}. Expressions are generated one node at a time along the pre-order traversal (from τ1 to τT ). For each node, a categorical distribution with parameters ψ defines the probabilities of selecting each node value from L. To capture the “context” of the expression as it is being generated, we condition this probability upon the selections of all previous nodes in that traversal. This conditional dependence can be achieved very generally using an RNN with parameters θ that outputs a probability vector ψ in autoregressive manner.
Specifically, the ith output vector ψ(i) of the RNN defines the probability distribution for selecting the ith node value τi, conditioned on the previously selected node values τ1:(i−1):
p(τi|τ1:(i−1); θ) = ψ (i) L(τi),
where L(τi) is the index in L corresponding to node value τi. The likelihood of the sampled expression is computed using the chain rule of conditional probability:
p(τ |θ) = |τ |∏ i=1 p(τi|τ1:(i−1); θ) = |τ |∏ i=1 ψ (i) L(τi)
The sampling process is illustrated in Figure 1 and described in Algorithm 1. Additional algorithmic details of the sampling process are described in Subroutines 1 and 2 in Appendix A. Starting at the root node, a node value is sampled according to ψ(1). Subsequent node values are sampled autoregressively in a depth-first, left-to-right manner until the tree is complete (i.e. all tree branches reach terminal nodes). The resulting sequence of node values is the tree’s pre-order traversal, which can be used to reconstruct the tree2 and its
1Given an expression tree (or equivalently, its pre-order traversal), the corresponding mathematical expression is unique; however, given an expression, its expression tree (or its corresponding traversal) is not unique. For example, x2 and x · x are equivalent expressions but yield different trees. For simplicity, we use τ somewhat abusively to refer to an expression where it technically refers to an expression tree (or equivalently, its corresponding traversal).
2In general, a pre-order traversal is insufficient to uniquely reconstruct the tree. However, in this context, we know how many child nodes each node has based on its value, e.g. “multiply” is a binary operator and thus has two children.
corresponding expression. Note that different samples of 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. Naively, the input to the RNN when sampling τi would be a representation (i.e. embedding or one-hot encoding) of the previously sampled token, τi−1. Indeed, this is typical in related autoregressive models, e.g. when generating sentences (Vaswani et al., 2017) or for neural architecture search (Zoph & Le, 2017). 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 in the expression tree. For example, the fifth and 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 node 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 a partial traversal is provided in Subroutine 2 in Appendix A.
Constraining the search space. Under our framework, it is straightforward to apply a priori constraints to 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 2 to prevent trivial expressions and a 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 a different constant. (3) The child of a unary operator should not be the inverse of that operator, e.g. log(exp(x)) is not allowed. (4) Direct descendants of trigonometric operators should not be trigonometric operators, e.g. sin(x + cos(x)) is not allowed because cosine is a 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 (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
For domains without this property, the number of children can be sampled from an additional RNN output. A pre-order traversal plus the corresponding number of children for each node is sufficient to uniquely reconstruct the tree.
Algorithm 1: Sampling an expression from the RNN 1 function SampleExpression(θ,L)
input : RNN with parameters θ; library of tokens L output: Pre-order traversal τ of an expression sampled from the RNN
2 τ = [] // Empty list 3 x = empty‖empty // Initial RNN input is empty parent and sibling 4 h0 = ~0 // Initialize RNN cell state to zero 5 for i = 1, . . . , T do 6 (ψ(i), hi) = RNN(x, hi−1; θ) 7 ψ(i) ← ApplyConstraints(ψ(i),L, τ) // Adjust probabilities 8 τi = Categorical(ψ(i)) // Sample next token 9 τ ← τ‖τi // Append token to traversal
10 if ExpressionComplete(τ) then 11 return τ 12 x← ParentSibling(τ) // Compute next parent and sibling 13 end 14 return τ
Algorithm 2: Deep symbolic regression 1 function DSR(α,N,L, X, y)
input : learning rate α; batch size N ; library of tokens L; input dataset (X, y) output: Best fitting expression τ?
2 Initialize RNN with parameters θ, defining distribution over expressions p(·|θ) 3 τ? = null 4 b = 0 5 repeat 6 T = {τ (i) ∼ p(·|θ)}i=1:N // Sample expressions (Algorithm 1) 7 T ← {OptimizeConstants(τ (i), X, y)}i=1:N // Optimize constants 8 R = {R(τ (i))− λCC(τ (i))}i=1:N // Compute rewards 9 ĝ = 1N ∑N i=1R(τ
(i))∇θ log p(τ (i)|θ) // Compute policy gradient 10 θ ← θ + α(ĝ1 + ĝ2) // Apply gradients 11 if maxR > R(τ?) then τ? ← τ (argmaxR) // Update best expression 12 return τ?
1 in Appendix A. This process ensures that all samples adhere to all constraints, without rejecting samples post hoc. In contrast, imposing constraints in GP-based symbolic regression can be problematic (Craenen et al., 2001). In practice, evolutionary operations that violate constraints are typically rejected post hoc (Fortin et al., 2012).

heading: 3.2 TRAINING THE RNN USING POLICY GRADIENTS
text: Optimizing the parameters of the sampled expressions. Once a pre-order traversal is sampled, we instantiate the corresponding symbolic expression. The expression may have several constant tokens, which can be viewed as model parameters. We train these model parameters by minimizing the mean-squared error
with respect to an input dataset using a nonlinear optimization algorithm, e.g. BFGS (Fletcher, 2013). We perform this inner optimization loop for each sampled expression before training the RNN.
Training the RNN using policy gradients. Given a distribution over mathematical expressions p(τ |θ) and a measure of performance of an expression R(τ), we consider the objective to maximize J(θ), defined as the expectation of R under expressions sampled from the distribution:
J(θ) ≡ Eτ∼p(τ |θ) [R(τ)]
We use REINFORCE (Williams, 1992) to maximize this expectation via gradient ascent:
∇θJ(θ) = ∇θEτ∼p(τ |θ) [R(τ)] = Eτ∼p(τ |θ) [R(τ)∇θ log p(τ |θ)]
This result allows us to estimate the expectation using samples from the distribution. Specifically, we can obtain an unbiased estimate of∇θJ(θ) by computing the sample mean over a batch of N sampled expressions T = {τ (i)}i=1:N :
∇θJ(θ) ≈ 1
N N∑ i=1 R(τ (i))∇θ log p(τ (i)|θ)
Reward function. A standard fitness measure in GP-based symbolic regression is normalized root-meansquare error (NRMSE), the root-mean-square error normalized by the standard deviation of the target values,
σy . That is, given a dataset of n number of (X, y) pairs, NRMSE = 1σy √ 1 n ∑n i=1(yi − ŷi)2, where ŷ = f(X) are the predicted values computed using the candidate expression f . Normalization by σy makes the metric commensurate across different datasets with potentially different ranges. However, metrics based on mean-square error exhibit extraordinarily large values for some expressions, e.g. an expression that incorrectly divides by an input variable with values near zero. For a gradient-based approach like DSR, this results in the gradient being dominated by the worst expressions, which can lead to instability. We found that a bounded reward function is more stable; thus, we applied a squashing function, yielding the reward function R(τ) = 1/(1 + NRMSE).3
We introduce the “vanilla” version of DSR in Algorithm 2. Below we describe several simple extensions.
Reward baseline. The above approximation to∇θJ(θ) is an unbiased gradient estimate, but in practice has high variance. To reduce variance, we include a baseline function b:
∇θJ(θ) ≈ 1
N N∑ i=1 [ R(τ (i))− b ] ∇θ log p(τ (i)|θ)
As long as the baseline is not a function of the current batch of expressions, the gradient estimate is still unbiased. We define the baseline function as an exponentially-weighted moving average of batches of rewards. Intuitively, the gradient step increases the likelihood of expressions above the baseline and decreases the likelihood of expressions below the baseline.
Complexity penalty. We include an optional complexity penalty that is added to the reward function. For simplicity, we consider the complexity metric |τ |, i.e. the number of nodes in the expression tree. More complicated metrics have been proposed that capture hierarchical features of the tree and/or deduced properties of the resulting expression (Vladislavleva et al., 2008).
3Since GP-based approaches using tournament selection only rely on the rankings of the fitness measure within the population, large fitness values are not problematic. SinceR(τ) is monotonic in NRMSE, GP is unaffected by squashing.
Algorithm 3: Deep symbolic regression with baseline, risk-seeking, entropy bonus, and complexity penalty
1 function DSR(α, β, λC , λH, ,N,L, X, y) input : learning rate α; moving average coefficient β; complexity coefficient λC ; entropy
coefficient λH; risk factor ; batch size N ; library of tokens L; input dataset (X, y) output: Best fitting expression τ?
2 Initialize RNN with parameters θ, defining distribution over expressions p(·|θ) 3 τ? = null 4 b = 0 5 repeat 6 T = {τ (i) ∼ p(·|θ)}i=1:N // Sample expressions (Algorithm 1) 7 T ← {OptimizeConstants(τ (i), X, y)}i=1:N // Optimize constants 8 R = {R(τ (i))− λCC(τ (i))}i=1:N // Compute rewards 9 R = (1− ) percentile ofR // Compute threshold
10 T ← {τ (i) : R(τ (i)) ≥ R } // Select subset of expressions 11 R ← {R(τ (i)) : R(τ (i)) ≥ R } // Select subset of rewards 12 ĝ1 = ReduceMean((R− b)∇θ log p(T |θ)) // Compute policy gradient 13 ĝ2 = ReduceMean(λH∇θH(T |θ)) // Compute entropy gradient 14 θ ← θ + α(ĝ1 + ĝ2) // Apply gradients 15 b← β · ReduceMean(R) + (1− β)b // Update baseline 16 if maxR > R(τ?) then τ? ← τ (argmaxR) // Update best expression 17 return τ?
Entropy bonus. We provide a bonus to the loss function proportional to the entropy of the sampled expressions. In accordance with the maximum entropy reinforcement learning framework (Haarnoja et al., 2018), this bonus serves two purposes. First, it encourages the RNN to explore more expressions, preventing premature convergence to a local optimum. In practice, this often leads to a better end result. Second, it encourages the RNN to assign equal likelihood to different expressions that have equal fitness.
Risk-seeking. The policy performance, J , is defined as an expectation. However, in practice, the performance of symbolic regression is measured by the single or few best expressions. Thus, we employ a novel risk-seeking technique in which only the top percentile samples from each batch are used in the gradient computation. This has the effect of increasing best-case performance at the expense of lower worst-case and average performances. This process is essentially the opposite of the EpOpt technique (Rajeswaran et al., 2016) used for risk-averse reinforcement learning, in which only the bottom percentile samples from each batch are used.
The complete algorithm, including reward baseline, complexity penalty, entropy bonus, and risk-seeking, is shown in Algorithm 3.

heading: 4 RESULTS AND DISCUSSION
text: Evaluating DSR. We evaluated DSR on a set of 12 commonly used symbolic regression benchmarks (Uy et al., 2011), as well as 4 additional variants in which we introduced real-valued constants to demonstrate the inner optimization loop. Each benchmark is defined by a ground truth expression, a training and testing dataset, and set of allowed operators, described in Table 2 in Appendix A. 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 at the end of training, and the ground truth function is used to determine whether the best found candidate expression was correctly recovered.
As a baseline, we compared against standard GP-based symbolic regression. To ensure fair comparison, the same constant optimizer (BFGS) was used for both methods. We ran independent training runs for GP and DSR for each benchmark expression (n = 100 for benchmarks without constants; n = 10 for benchmarks with constants). For each experiment, we generated 1,000 candidate expressions per generation/iteration for 1,000 generations/iterations, resulting in 1,000,000 total expressions. For each training run, the expression with the best reward is selected and we record the NRMSE on the test data.
For GP, we used the open-source software package “deap” (Fortin et al., 2012). For DSR, the RNN comprised a single-layer LSTM of 32 hidden units. Additional hyperparameters and experiment details are provided in Appendix A.
In Table 1, we report the percentage of runs that correctly recover the expression and NRMSE on the test data for each benchmark. DSR significantly outperforms GP in its ability to exactly recover benchmark expressions. DSR also outperforms GP in the average NRMSE across all expressions; however, we observe that for the few expressions with low or zero recovery rate (e.g. Nguyen-7, Nguyen-8, and Nguyen-12), GP sometimes exhibits lower NRMSE. One explanation is that GP is more prone to overfitting the expression to the dataset. As an 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 scratch each iteration after learning from a gradient update, making it less prone to overfitting.
Surprisingly, DSR consistently performed best without a complexity penalty, i.e. λC = 0. Due to the autoregressive nature of the RNN, shorter expressions tend to exhibit higher likelihood than longer ones. We postulate that this property produces a self-regularization effect that precludes the need for an explicit complexity penalty.
Ablation studies. Algorithm 3 includes several additional components relative to the “vanilla” Algorithm 2. We performed a series of ablation studies to quantify the effect of each of these components, along with the effects of the various constraints on the search space, and including the parent and sibling as input to the RNN instead of the previous node value. In Figure 2, we performed DSR on the set of 12 Nguyen benchmarks for each ablation. DSR is still competitive with GP even when removing all improvements and all constraints.
Noisy data and amount of data. We evaluated the robustness of DSR to noisy data by adding independent Gaussian noise to the dependent variable, with mean zero and standard deviation proportional to the root-mean-square of the dependent variable in the training data. In Figure 3, we varied the proportionality constant from 0 (noiseless) to 10−1 and compared the performance of GP and DSR across the set of 12 Nguyen benchmarks. DSR still outperforms GP in both recovery rate and NRMSE across noise levels.
Symbolic regression excels in the low-data setting when data is noiseless, hence, the benchmark expressions included herein include only 20 data points (see Table 2). With added noise, increasing the amount of data
smooths the reward function and may help prevent overfitting. Thus, we repeated the noise experiments using the same benchmarks but with 10-fold larger training datasets (200 points data points). As expected, recovery rates tend to increase for both methods; however, DSR maintains a much larger improvement than GP at higher noise levels.

heading: 5 CONCLUSION
text: We introduce an unconventional approach to symbolic regression based on reinforcement learning that outperforms a standard GP-based method on recovering exact expressions on benchmark problems, both with and without added noise. Since both DSR and GP generate expression trees, there are many opportunities for hybrid methods, for example including several generations of evolutionary operations in the inner optimization loop. From the perspective of AutoML, the main contributions are defining a flexible distribution over hierarchical, variable-length objects that allows imposing in situ constraints, and using risk-seeking training to optimize best-case performance. Thus, we note that our framework is easily extensible to domains outside symbolic regression, which we save for future work; for example, searching the space of organic molecular structures for high binding affinity to a reference compound. We chose symbolic regression to demonstrate our framework in part because of the large search space, broad applicability, computationally expedient inner optimization loop (sub-second), and availability of vetted benchmark problems and baseline methods.

heading: APPENDIX A
text: Hyperparameters. DSR hyperparameters are listed in Table 3. GP hyperparameters are listed in Table 4. The same hyperparameters were used for all experiments and all benchmark expressions.
Additional details for performance comparison experiments. Details of the benchmark symbolic regression problems are shown in Table 2. All benchmarks use the function set {+,−,×,÷, sin, cos, exp, log}. To ensure closures, we use protected versions of operators: log returns the logarithm of the absolute value of its argument, and ÷, exp, and log return 1 for arguments that would cause overflow or other errors. Benchmarks without constants can be recovered exactly, thus recovery is defined by exact correctness (modulo floating point precision error). Note Nguyen-8 can be recovered via exp( xx+x log(x)) and Nguyen-11 can be recovered via exp(y log(x)).
Table 3: DSR hyperparameters
Parameter Value Batch size 1,000 Iterations 1,000
Learning rate (α) 0.0003 Entropy coefficient (λH) 0.08
Complexity coefficient (λC) 0 Moving average coefficient (β) 0.5
Risk factor ( ) 0.1
Table 4: GP hyperparameters
Parameter Value Population size 1,000
Generations 1,000 Fitness function NRMSE
Initialization method Full Selection type Tournament Tournament size (k) 3 Crossover probability 0.5 Mutation probability 0.1
Minimum subtree depth (dmin) 0 Maximum subtree depth (dmax) 2
For benchmarks with constants, constants are optimized using BFGS with an initial guess of 1.0 for each constant. 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 floating point constants cannot be recovered exactly, for benchmarks with constants we manually determined correctness of the functional form by inspection. Since constant optimization is a computational bottleneck, we limited each expression to three constants for both DSR and GP experiments.
For GP, 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.
Additional details for ablation studies. In Figure 2, “Parent/sibling” denotes that the previous node of the traversal is provided as input to the RNN, rather than the parent and sibling node. “Risk-seeking” denotes no risk-seeking, equivalent to = 1. “Entropy bonus” denotes no entropy bonus, equivalent to λH = 0. “Reward baseline” denotes no reward baseline, equivalent to β = 0. “All improvements” denotes combining ablations for Parent/sibling, Risk, Entropy, and Baseline. “Constrain trig” denotes no constraint precluding nested trigonometric operators. “Constrain inverse” denotes no constraint precluding inverse unary operators. “Constrain min/max” denotes no constraint precluding minimum or maximum length. (If the maximum length is reached, the expression is appended with x until complete.) “Constraints” denotes combining ablations for Trig, Inverse, and Min/max. “All constraints & improvements” denotes combining all ablations.
Training curves. Figures 4 and 5 show the reward (1/(1 + NRMSE)) and recovery rate, respectively, as a function of training step (DSR iteration or GP generation). For benchmarks with constants, the constant optimizer can allow both algorithms to quickly reach reward near 1.0. For these benchmarks, we provide zoomed inset plots demonstrate if and when all independent training runs correctly recovery the expression (Figure 4: Constant-1, Constant-2, and Constant-3). Note that the NRMSE and recovery values in Table 1 correspond to the final point on each curve in Figures 4 and 5, respectively.
Additional subroutines. DSR includes several subroutines used when sampling an expression from the RNN and during training. In Subroutine 1, we describe the function ApplyConstraints(ψ,L, τ) used in Algorithm 1, which zeros out the probabilities of tokens that would violate any given constraints. Within this subroutine, the user-specific function ViolatesConstraint(τ,Li) returns TRUE if adding the ith token of the library L to the partial traversal τ would violate any user-specified constraints, and FALSE otherwise. In Subroutine 2, we describe the function ParentSibling(τ) used in Algorithm 1, which computes the parent and sibling of the next token to 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 is the parent and there is no sibling. Otherwise, the subroutine iterates backward through the traversal until finding a node with an unselected child node. That node is the parent and the subsequent node is the child. Within this subroutine, the function Arity(τi) simply returns the arity (number of arguments) of token τi, i.e. two for binary operators, one for unary operators, or zero for input variables and constants.
In Subroutine 3, we describe the function OptimizeConstants(τ,X, y) used in Algorithms 2 and 3, which optimizes the placeholder constants c of an expression with respect to input dataset (X, y) using a black-box optimizer, e.g. BFGS. Within this subroutine, the function Instantiate(τ) instantiates the symbolic expression as a function f(X; c) with inputs X and parameters (constants) c, and the function ReplaceConstants(τ, c?) replaces the placeholder constants in the expression with the optimized constants.
Subroutine 1: Apply generic constraints in situ when sampling from the RNN 1 function ApplyConstraints(ψ,L, τ)
input : Categorical probabilities ψ; corresponding library of tokens L; partially sampled traversal τ output: Adjusted categorical probabilities ψ
2 L = |ψ| // Length of library 3 for i = 1, . . . , L do 4 if ViolatesConstraint(τ,Li) then 5 ψi ← 0 // Constrain that token 6 end 7 ψ ← ψ∑
i ψi // Normalize back to 1
8 return ψ
Subroutine 2: Computing parent and sibling inputs to the RNN 1 function ParentSibling(τ)
input : Partially sampled traversal τ output: Parent and sibling tokens of the next token to be sampled
2 L = |τ | // Length of partial traversal 3 counter = 0 // Counter for number of unselected nodes 4 if Arity(τL) > 0 then 5 parent = τL 6 sibling = empty 7 return parent‖sibling 8 for i = L, . . . , 1 do // Iterate backward 9 counter← counter + Arity(τi)− 1
10 if counter = 0 then 11 parent = τi 12 sibling = τi+1 13 return parent‖sibling 14 end
Subroutine 3: Optimize the constants of an expression (inner optimization loop) 1 function OptimizeConstants(τ,X, y)
input : Expression τ with placeholder constants c; input dataset (X, y) output: Expression τ? with optimized constants c?
2 f(X, c) = Instantiate(τ) // Instantiate the symbolic expression 3 c? = argminc ‖y − f(X, c)‖22 // Minimize error (e.g. with BFGS) 4 τ? = ReplaceConstants(τ, c?) // Replace placeholder constants 5 return τ?

Review summaries:
This paper presents deep symbolic regression (DSR), which uses a recurrent neural network to learn a distribution over mathematical expressions and uses policy gradient to train the RNN for generating desired expressions given a set of points. The RNN model is used to sample expressions from the learned distribution, which are then instantiated into corresponding trees and evaluated on a dataset. The fitness on the dataset is used as the reward to train the RNN using policy gradient. In comparison to GP, the presented DSR approach recovers exact symbolic expressions in majority of the benchmarks.

This paper presents a RNN-RL based method for the symbolic regression problem. The problem is new (to Deep RL) and interesting. My main concern is about the proposed method, where the three RL related equations (not numbered) at page 5 are also direct copy-from-textbook policy gradient equations without specific adaptation to the new application considered in this paper, which is very strange. The two conditional probability definitions considered at page 3 are not mentioned in later text. These are only fractions of the underlying method and by reading the paper back and forth several times, it is not clear of the basic algorithmic flowchart, let alone more detailed description of the related parameters. Without these information, it is impossible to have a fair judge of the novelty and feasibility of the proposed method. The empirical results are also limited in small dataset, which makes it hard to verify the generality of the superior claim.