Title: NEURAL-GUIDED SYMBOLIC REGRESSION WITH ASYMPTOTIC CONSTRAINTS

Abstract: Symbolic regression is a type of discrete optimization problem that involves searching expressions that fit given data points. In many cases, other mathematical constraints about the unknown expression not only provide more information beyond just values at some inputs, but also effectively constrain the search space. We identify the asymptotic constraints of leading polynomial powers as the function approaches 0 and ∞ as useful constraints and create a system to use them for symbolic regression. The first part of the system is a conditional production rule generating neural network which preferentially generates production rules to construct expressions with the desired leading powers, producing novel expressions outside the training domain. The second part, which we call Neural-Guided Monte Carlo Tree Search, uses the network during a search to find an expression that conforms to a set of data points and desired leading powers. Lastly, we provide an extensive experimental validation on thousands of target expressions showing the efficacy of our system compared to exiting methods for finding unknown functions outside of the training set.

1. INTRODUCTION:

The long standing problem of symbolic regression tries to search expressions in large space that fit given data points (Koza & Koza, 1992; Schmidt & Lipson, 2009) . These mathematical expressions are much more like discovered mathematical laws that have been an essential part of the natural sciences for centuries. Since the size of the search space increases exponentially with the length of expressions, current search methods can only scale to find expressions of limited length. Moreover, current symbolic regression techniques fail to exploit a key value of mathematical expressions that has traditionally been well used by natural scientists. Symbolically derivable properties such as bounds, limits, and derivatives can provide significant guidance to finding an appropriate expression.

In this work, we consider one such property corresponding to the behavior of the unknown function as it approaches 0 and ∞. Many expressions have a defined leading polynomial power in these limits. For examples when x → ∞, 2x 2 + 5x has a leading power of 2 (because the expression behaves like x 2 ) and 1/x 2 + 1/x has a leading power of -1. We call these properties "asymptotic constraints" because this kind of property is known a priori for some physical systems before the detailed law is derived. For example, most materials have a heat capacity proportional to T 3 at low temperatures T and the gravitational field of planets (at distance r) should behave as 1/r as r → ∞.

Asymptotic constraints not only provide more information about the expression, leading to better extrapolation, but also constrain the search in the desired semantic subspace, making the search more tractable in much larger space. These constraints can not be simply incorporated using syntactic restrictions over the grammar of expressions. We present a system to effectively use asymptotic constraints for symbolic regression, which has two main parts. The first is a conditional production rule generating neural network (NN) of the desired polynomial leading powers that generates production rules to construct novel expressions (both syntactically and semantically) and, more surprisingly, generalize to leading powers not in the training set. The second part is a Neural-Guided Monte Carlo Tree Search (NG-MCTS) that uses this NN to probabilistically guide the search at every step to find expressions that fit a set of data points.

Finally, we provide an extensive empirical evaluation of the system compared to several strong baseline techniques. We examine both the NG-MCTS and conditional production rule generating NN alone. In sharp contrast to almost all previous symbolic regression work, we evaluate our technique on thousands of target expressions and show that NG-MCTS can successfully find the target expressions in a much larger fraction of cases (71%) than other methods (23%) with search space sizes of more than 10 50 expressions.

In summary, this paper makes the following key contributions: 1) We identify asymptotic constraints as important additional information for symbolic regression tasks. 2) We develop a conditional production rule generating NN to learn a distribution over (syntactically-valid) production rules conditioned on the asymptotic constraints. 3) We develop the NG-MCTS algorithm that uses the conditional production rule generating NN to efficiently guide the MCTS in large space. 4) We extensively evaluate our production rule generating NN to demonstrate generalization for leading powers, and show that the NG-MCTS algorithm significantly outperforms previous techniques on thousands of tasks.

2. PROBLEM DEFINITION:

In order to demonstrate how prior knowledge can be incorporated into symbolic regression, we construct a symbolic space using a context-free grammar G:

O → S S → S'+'T | S'-'T | S'*'T | S'/'T | T T → '('S')' | 'x' | '1' (1)

This expression space covers a rich family of rational expressions, and the size of the space can be further parameterized by a bound on the maximum sequence length. For an expression f (x), the leading power at x 0 is defined as P x→x0 [f ] = p s.t. lim x→x0 f (x)/x p = non-zero constant. In this paper, we consider the leading powers at x 0 ∈ {0, ∞} as additional specification.

Let S(G, k) denote the space of all expressions in the Grammar G with a maximum sequence length k. Conventional symbolic regression searches for a desired expression f (x) in the space of expressions S(G, k) that conforms to a set of data points {(x, f (x)) | x ∈ D train }, i.e. find a g(x) ∈ S(G, k) : φ(g(x), D train ), where φ denotes the acceptance criterion, usually root mean square error (RMSE) . With the additional specification of leading powers c (0) and c (∞) at 0 and ∞, the problem becomes: find a g(x) ∈ S(G, k) :

φ(g(x), D train )∧(P x→0 [g] = c (0) )∧(P x→∞ [g] = c (∞) ).

3. CONDITIONAL PRODUCTION RULE GENERATING NEURAL NETWORK:

It is difficult to directly incorporate the asymptotic constraints as syntactic restrictions over the grammar. We, therefore, develop a neural architecture that learns to generate production rule in the grammar conditioned on the given asymptotic constraints.

GRU:

⨉ + + + + O → S S → S '/' T S → T T → '1' O → S S → S '+' T S → S '-' T S → S '*' T S → S '/' T S → T T → '(' S ')' T → 'x' T → '1'

leading power x→0 leading power x→∞ Figure 1 (a) and (c) show an example of how an expression is parsed as a parse tree by the grammar defined in Eq. ( 1). The parse tree in Figure 1 (c) can be serialized into a production rule sequence r 1 , . . . , r L by a preorder traversal (Figure 1 (d)), where L denotes the length of the production rule sequence. Figure 1 (b) shows the leading powers of the exemplary expression in Figure 1 (a). The conditional distribution of an expression is parameterized as a sequential model

EQUATION

EQUATION

We build a NN (as shown in Figure 1 (e)) to predict the next production rule r t+1 from a partial sequence r 1 , . . . , r t and conditions c (0) , c (∞) . During training, each expression in the training set is first parsed as a production rule sequence. Then a partial sequence of length t ∈ {1, . . . , L -1} is sampled randomly as the input and the (t + 1)-th production rule is selected as the output (see blue and orange text in Figure 1(d) ). Each production rule of the partial sequence is represented as an embedding vector of size 10. The conditions are concatenated with each embedding vector. This sequence of embedding vectors are fed into a bidirectional Gated Recurrent Units (GRU) (Cho et al., 2014) with 1000 units. A softmax layer is applied to the final output of GRU to obtain the raw probability distribution over the next production rules in Eq. ( 1).

Note that not all the production rules are grammatically valid as the next production rule for a given partial sequence. The partial sequence is equivalent to a partial parse tree. The next production rule expands the leftmost non-terminal symbol in the partial parse tree. For the partial sequence colored in blue in Figure 1 (d), the next production rule expands non-terminal symbol T , which constrains the next production rule to only those with left-hand-side symbol T . We use a stack to keep track of non-terminal symbols in a partial sequence as described in GVAE (Kusner et al., 2017) . A mask of valid production rules is computed from the input partial sequence. This mask is applied to the raw probability distribution and the result is normalized to 1 as the output probability distribution.

The training loss is calculated as the cross entropy between the output probability distribution and the next target production rule. It is trained from partial sequences sampled from expressions in the training set using validation loss for early stopping.1 

4. NEURAL-GUIDED MONTE CARLO TREE SEARCH:

We now briefly describe the NG-MCTS algorithm that uses the conditional production rule generating NN to guide the symbolic regression search. The discrepancy between the best found expression g(x) and the desired f (x) is evaluated on data points and leading powers. The error on data points is measured by RMSE ∆g 

{•} = x∈D {•} (f (x) -g(x)) 2 /|D

] = |P x→0 [f ] -P x→0 [g]| + |P x→∞ [f ] -P x→∞ [g]|.

The default choice of objective function for symbolic regression algorithms is ∆g train alone. With additional leading powers constraint, the objective function can be defined as ∆g train + ∆P [g], which minimizes both the RMSE on the training points and the absolute difference of the leading powers.

Most symbolic regression algorithms are based on EA (Schmidt & Lipson, 2009) , where it is nontrivial to incorporate our conditional production rule generating NN to guide the generation strategy in a step-by-step manner, as the mutation and cross-over operators perform transformations on fully completed expressions. However, it is possible to incorporate a probability distribution over expressions in many heuristic search algorithms such as Monte Carlo Tree Search (MCTS). MCTS is a heuristic search algorithm that has been shown to perform exceedingly well in problems with large combinatorial space, such as mastering the game of Go (Silver et al., 2016) and planning chemical syntheses (Segler et al., 2018) . In MCTS for symbolic regression, a partial parse tree sequence r 1 , . . . , r t can be defined as a state s t and the next production rule is a set of actions {a}. In the selection step, we use a variant of the PUCT algorithm (Silver et al., 2016; Rosin, 2011) for exploration. For MCTS, the prior probability distribution p(a i |s t ) is uniform among all valid actions.

We develop NG-MCTS by incorporating the conditional production rule generating NN into MCTS for symbolic regression. Figure 1 (f) presents a visual overview of NG-MCTS. In particular, the prior probability distribution p(a i |s t , c (0) , c (∞) ) is computed by our conditional production rule generating NN on the partial sequence and the desired conditions. We run MCTS for 500 simulations for each desired expression f (x). The exploration strength is set to 50 and the production rule sequence length limit is set to 100. The total number of expressions in this combinatorial space is 3 × 10 93 . We run evolutionary algorithm (EA) (Appendix G) and grammar variational autoencoder (GVAE) (Kusner et al., 2017) (Appendix H) with comparable computational setup for comparison.

5. EVALUATION:

Dataset and Search Space: We denote the leading power constraint as a pair of integers

(P x→0 [f ], P x→∞ [f ]

) and define the complexity of a condition by

M [f ] = |P x→0 [f ]| + |P x→∞ [f ]|. Obviously, expressions with M [f ] = 4 are more complicated to construct than those with M [f ] = 0.

We create a dataset balanced on each condition, as described in Appendix A. The conditional production rule generating NN is trained on 28837 expressions and validated on 4095 expressions, both with 

M [f ] ≤ 4.

5.1. EVALUATION OF SYMBOLIC REGRESSION TASKS:

We now present the evaluation of our NG-MCTS method, where each step in the search is guided by the conditional production rule generating NN. Recent developments of symbolic regression methods (Kusner et al., 2017; Sahoo et al., 2018) compare methods on only a few expressions, which may cause the performance to depend on random seed for initialization and delicate tuning. To mitigate this, we apply different symbolic regression methods to search for expressions in holdout sets with thousands of expressions and compare their results in Table 1 . Additional comparison on a subset of holdout sets is reported in Appendix H. In order to examine the effect of leading powers, we use leading powers alone in MCTS (PW-ONLY).

The median of ∆P [g] for hard expressions is reduced to 1 but the medians of ∆g int. and ∆g ext. are significantly higher. We then add leading powers to the objective function together with data points. MCTS + PW does not have a notable difference to MCTS. However, EA + PW improves solved expressions to 23.37% and ∆P [g] of hard expressions is 0. This indicates adding leading power constraints in the objective function is helpful for symbolic regression. Most importantly, we observe step-wise guidance of NN conditioned on leading powers can lead to even more significant improvements compared to adding them in the objective function. NG-MCTS solves 71.22% expressions in the holdout set, three times over the best EA + PW. Note that both MCTS + PW and EA + PW have access to the same input information. Although EA has the lowest medians of ∆g train and ∆g int. , NG-MCTS is only slightly worse. On the other hand, NG-MCTS outperforms on ∆g ext.

and ∆P [g], which indicates that step-wise guidance of leading powers helps to generalize better in extrapolation than all other methods.

We also apply the aforementioned methods to search expressions in holdout sets 

M [f ] = 5 and M [f ] = 6.

NG-MCTS 1 -x + (1/x) + x × x 0.00 0.00 0.00 0 GVAE (x) -(1/x)/(x × x/x) + x 0.47 0.29 34.9 2 GVAE + PW ((1/x) -x + x) -((1 -x) × x) 1.0 1.0 1.0 0 EA (x + x) 0.52 0.46 34.8 3 EA + PW ((1/x) + (x × x))

1.15 1.10 6.16 0 of asymptotic constraints. We propose a force field potential U (x) = 1/x + x + (x -1) 2 with Coulomb interaction, uniform electric field and harmonic interaction. Assuming the true potential is unknown, the goal is to discover this expression. As a physical potential, besides values at D train : {1.2, 1.6, 2.0, 2.4, 2.8}, researchers also know the short (x → 0) and long range (x → ∞) behaviors as leading powers. Table 2 shows the expressions found by NG-MCTS, GVAE, GVAE + PW, EA and EA + PW, which are plotted in Figure 3 . NG-MCTS can find the desired expression.

The second best method is GVAE + PW, differing by a constant of 1 from the true expression.

5.3. EVALUATION OF CONDITIONAL PRODUCTION RULE GENERATING NN:

NG-MCTS significantly outperforms other methods on searching expressions in large space. To better understand the effective guidance from NN, we demonstrate its ability to generate syntactically and semantically novel expressions given desired conditions. In order to examine the NN alone, we directly sample from the model by Eq. ( 2) instead of using MCTS. The model predicts the probability distribution over the next production rules from the starting rule r 1 : O → S and desired condition (c (0) , c (∞) ). The next production rule is sampled from distribution p θ (r 2 |r 1 , c (0) , c (∞) ) and then appended to r 1 . Then r 3 is sampled from p θ (r 3 |r 1 , r 2 , c (0) , c (∞) ) and appended to [r 1 , r 2 ]. This procedure is repeated until [r 1 , . . . , r L ] form a parse tree where all the leaf nodes are terminal, or the length of generated sequence reaches the prespecified limit, which is set to 100 for our experiments.

Baseline Models We compare NN with a number of baseline models that provide a probability distribution over the next production rules. All these distributions are masked by the valid production rules computed from the partial sequence before sampling. For each desired condition within

|P x→0 [f ]| ≤ 9 and |P x→∞ [f ]| ≤ 9, k = 100 expressions are generated.

We consider the following baseline models: i) Neural Network No Condition (NNNC): same setup as NN except no conditioning on leading powers, ii) Random: uniform distribution over valid next production rules, iii) Full History (FH): using full partial sequence and conditions to sample next production rule from its empirical distribution, iv) Full History No Condition (FHNC), v) Limited History (LH) (l): sampling next production rule from its empirical distribution given only the last l rules in the partial sequence, and vi) Limited History No Condition (LHNC) (l). Note that the aforementioned empirical distributions are derived from the training set {f }. For limited history models, if l exceeds the length of the partial sequence, we instead take the full partial sequence. More details about the baseline models can be found in Appendix B.

Metrics We propose four metrics to evaluate the performance. For each condition (c ∞) ). i) Success Rate: proportion of generated expressions with leading powers

EQUATION

(P x→0 [g i ], P x→∞ [g i ]) that match the desired condition, ii) Mean L1-distance: the mean L1-distance between (P x→0 [g i ], P x→∞ [g i ]) and (c (0) , c (∞) ), iii) Syntactic

Novelty Rate: proportion of generated expressions that satisfy the condition and are syntactically novel (no syntactic duplicate of generated expression in the training set), and iv) Semantic Novelty Rate: proportion of expressions satisfying the condition and are semantically novel (the expression obtained after normalizing the generated expression does not exist in the training set). For example, expressions x + 1 and (1) + x are syntactically different, but semantically duplicate. We use simplify function in SymPy (Meurer et al., 2017) to normalize expressions. To avoid inflating the rates of syntactic and semantic novelty, we only count the number of unique syntactic and semantic novelties in terms of their expressions and simplified expressions, respectively. Model Furthermore, the L1-distance measures the deviation from the desired condition when not matching exactly. NN has the least mean L1-distance in the out-of-sample condition region. This suggests that for the unmatched expressions, NN prefers expressions with leading powers closer to the desired condition than all other models. NN outperforms the other models not only on the metrics aggregated over all conditions, but also for individual conditions. Figure 4 shows the metrics for NN and LHNC (8) on each condition. NN performs better in the in-sample region (inside the red boundary) and also generalizes to more conditions in the out-of-sample region (outside the red boundary). Qualitative Evaluation To better comprehend the learned NN model and its generative behavior, we also perform a task of expression completion given a structure template of the form 1/ and a variety of desired conditions in Table 4 . For each condition, 1000 expressions are generated by NN and the probability of each syntactically unique expression is computed from its occurrence. We first start with c (0) = 0, c (∞) = 1. The completed expression g(x) is required to be a nonzero constant as x → 0 and g(x) → x as x → ∞. The top three probabilities are close since this task is relatively easy to complete. We then repeat the task on c (0) = -1, c (∞) = 1 and c (0) = -2, c (∞) = 2, which are still in the in-sample condition region and the model can still complete expressions that match the desired conditions. We also show examples of c (0) = -3, c (∞) = 2, which is in the out-of-sample condition region. Note that to match condition c (∞) = -3, more complicated completion such as 1/(x * (x * x)) has to be constructed by the model. Even in this challenging case, the model can still generate some expressions that match the desired condition. We also show examples of the syntactic novelties learned by model in Appendix L. The EA varies the structures to search new expressions using a score that accounts for both accuracy and the complexity of the expression. This approach has been further used to get empirical expressions in electronic engineering (Ceperic et al., 2014) , water resources (Klotz et al., 2017) , and social science (Truscott & Korns, 2014) . GVAE (Kusner et al., 2017) was recently proposed to learn a generative model of structured arithmetic expressions and molecules, where the latent representation captures the underlying structure. This model was further shown to improve a Bayesian optimization based method for symbolic regression.

Quantitative Evaluation:

M [f ] ≤ 4 M [f ] > 4 M [f ] ≤ 4 = 5 = 6 =7 Syn (%)

6. RELATED WORK:

M [f ] c (0) c (∞) Expression Probability Match 1 0 1 1/ (1 + 1) -x 7.8% √ 1/ (1 + x) -x 7.7% √ 1/ x -x 7.0% × 2 -1 1 1/ (x + x) -x 17.3% √ 1/ (x) -x 12.0% √ 1/ x -x 6.8% √ 4 -2 2 1/ (x * x) -(x * x) 48.5% √ 1/ x -(x * x) 12.5% × 1/ (x * x) -x 7.6% × 5 -3 2 1/ (x * x) -(x * x) 29.5% × 1/ (x * x) -(x * (x * x)) 19.3% × 1/ (x * (x * x)) -(x *

Similar to these approaches, most other approaches search for expressions from scratch using only data points (Schmidt & Lipson, 2009; Ramachandran et al., 2017; Ouyang et al., 2018) without other symbolic constraints about the desired expression. Abu-Mostafa (1994) suggests incorporating prior knowledge of a similar form to our property constraints, but actually implements those priors by adding additional data points and terms in the loss function.

Neural Program Synthesis: Program synthesis is the task of learning programs in a domain-specific language (DSL) that satisfy a given specification (Gulwani et al., 2017) . It is closely related to symbolic regression, where the DSL can be considered as a grammar defining the space of expressions. Some recent works use neural networks for learning programs (Devlin et al., 2017; Balog et al., 2016; Parisotto et al., 2017; Vijayakumar et al., 2018) . RobustFill (Devlin et al., 2017) trains an encoderdecoder model that learns to decode programs as a sequence of tokens given a set of input-output examples. For more complex DSL grammars such as Karel that consists of nested control-flow, an additional grammar mask is used to ensure syntactic validity of the decoded programs (Bunel et al., 2018) . However, these approaches only use examples as specification. Our technique can be useful to guide search for programs where the program space is defined using grammars such as SyGuS with additional semantic constraints other than only examples (Alur et al., 2013; Si et al., 2018) .

7. CONCLUSION:

We present a two-step framework of first learning a neural network of the relation between syntactic structure and leading powers and then using that to guide MCTS for efficient search. This framework is evaluated in the context of symbolic regression and focused on the leading power properties on thousands of desired expressions, a much larger evaluation set than benchmarks considered in existing literature. We plan to further extend the applicability of this framework to cover other symbolically derivable properties of expressions. Similar modeling ideas could be equally useful in general program synthesis settings, where other properties such as the desired time complexity or maximum control flow nesting could be used as constraints.

A DATASET:

To create an initial dataset, we first enumerate all possible parse trees from

O → S S → S'+'T | S'-'T | S'*'T | S'/'T | T T → '('S')' | 'x' | '1'

within ten production rules. Then we repeat the following downsampling and augmentation operations for four times to expand the dataset for longer expressions and diverse conditions.

Downsampling Expressions are grouped by their simplified expressions computed by SymPy (Meurer et al., 2017) . We select 20 shortest expressions in terms of string length from each group. These expressions are kept to ensure syntactical diversity and avoid having too long expressions.

Augmentation For each kept expression, five new expressions are created by randomly replacing one of the 1 or x symbols by

(1/x), (x/(1 + x)), (x/(1 -x)), (1/(1 + x)), (1/(1 -x)), (1 -x), (1 + x), (x * x), (x * (1 + x)), (x * (1 -x))

. These five newly created expressions are added back to the pool of kept expressions to form an expanded set of expressions.

After repeating the above operations for four times, we apply the downsampling step again in the end.

To make the dataset balanced on each condition, we keep 1000 shortest expressions in terms of string length for each condition. In this way, we efficiently create an expanded dataset which is not only balanced on each condition but also contains a large variety of expressions with much longer length than the initial dataset. Compared to enumerating all possible expressions given the maximum length, the created dataset is much sparser and smaller. A .1 shows the minimum, median and maximum of the lengths of the production rule sequence. The last column, space size, shows the number of possible expressions within the maximum length of the production rule sequence (including those nonterminal expressions). It is computed recursively as follows. Let N * ,i denote the number of possible expressions with length ≤ i and whose production rule sequences start with symbol * , where * can be O, S and T . Then we have

N S,i = 4 i-1 p=0 (N S,p • N T,i-1-p ) + N T,i-1 ;

N T,i = N S,i-1 + 2, and the initial condition is N S,0 = N T,0 = 1. We can obtain N S,i using the recursive formula. The quantity of interest N O,i = N S,i-1 given that the first production rule is pre-defined as O → S. Holdout sets M [f ] = 5, 6 not only contain expressions of higher leading powers but also of longer length, which is challenging for generalization both semantically and syntactically. 

B DETAILS OF THE BASELINE MODELS:

We provide more details of the baseline models we proposed to be compared with our NN model. Using the same notation as in Eq. ( 2), the conditional distribution of the next production rule given the partial sequence and the desired condition is denoted by p(r t+1 |r 1 , . . . , r t , c (0) , c (∞) ).

The baseline models are essentially different ways to approximate the conditional distribution using empirical distributions.

Limited History (LH) (l):

The conditional distribution is approximated by the empirical conditional distribution given at most the last l production rules of the partial sequence and the desired condition. We derive the empirical conditional distribution from the training set by first finding all the partial sequences therein that match the given partial sequence and desired condition. Then we compute the proportion of each production rule that appears as the next production rule of the found partial sequences. The proportional is therefore the empirical conditional distribution. To avoid introducing an invalid next production rule, the empirical conditional distribution is multiplied by the production rule mask of valid next production rules, and renormalized.

p LH(l) = p(r t+1 |r t-l+1 , . . . , r t , c (0) , c (∞) ).

Full History (FH) The conditional distribution is approximated by the empirical conditional distribution given the full partial sequence and the desired condition. The empirical conditional distribution is derived from the training set similarly as the LH model.

p FH = p(r t+1 |r 1 , . . . , r t , c (0) , c (∞) ).

Limited History No Condition (LHNC) (l):

The conditional distribution is approximated by the empirical conditional distribution given at most the last l production rules of the partial sequence only, where the desired condition is ignored. The empirical conditional distribution is derived from the training set similarly as the LH model.

p LHNC(l) = p(r t+1 |r t-l+1 , . . . , r t ).

Full History No Condition (FHNC):

The conditional distribution is approximated by the empirical conditional distribution given the full partial sequence only, where the desired condition is ignored. The empirical conditional distribution is derived from the training set similarly as the LH model. p FHNC = p(r t+1 |r 1 , . . . , r t ).

C DETAILS OF AGGREGATING METRICS ON DIFFERENT LEVELS:

The metrics are aggregated on different levels. We compute the average of success rates, syntactic novelty rates and semantic novelty rates over all the conditions in the in-sample condition region. The out-of-sample condition region is not bounded, and hence we consider the region within

|P x→0 [f ]| ≤ 9 and |P x→∞ [f ]| ≤ 9.

Since the average can be arbitrarily small if the boundary is arbitrarily large, instead we compute the total number of success, syntactic novelty and semantic novelty expressions over all the conditions in the out-of-sample region.

The L1-distance is not well defined for an expression which is non-terminal or ∞. For both cases, we specify the L1-distance as 18, which is the L1-distance between conditions (0, 0) and (9, 9). A model that generates expressions with a high success rate (i.e., satisfying the desired condition most of the times) but lacking of diversity is problematic. To demonstrate the diversity of expressions generated by our NN model, we generate 1000 expressions on each desired condition, and compute the cumulative counts of unique expressions and unique simplified expressions that satisfy the desired condition among the first number of expressions of the 1000 generated expressions, respectively. Figure D.1 shows the cumulative counts on three typical desired conditions (c (0) , c (∞) ) = (0, 0), (5, 5), (-2, -3). We can observe that the counts steadily increase as the number of generated expressions under consideration increases. Even with 1000 expressions, the counts have not been saturated.

D DIVERSITY OF GENERATED EXPRESSIONS:



E ADDITIONAL PLOTS OF METRICS FOR PRODUCTION RULE GENERATING NN AND BASELINE MODELS:

Due to the space limit of the paper, Figure 4 in the main text only contains the plots of metrics on each condition for LHNC (8) and NN. Additional plots of all the baseline models in 

F SELECTION STEP IN MONTE CARLO TREE SEARCH:

In the selection step, we use a variant of the PUCT algorithm (Silver et al., 2016; Rosin, 2011) for exploration,

EQUATION

where N (s t , a) is the number of visits to the current node, b N (s t , b) is the number of visits to the parent of the current node, c puct controls the strength of exploration and P (s t , a) is the prior probability of action a. This strategy initially prefers an action with high prior probability and low visits for similar tree node quality Q(s t , a).

G EVOLUTIONARY ALGORITHM:

We implemented the conventional symbolic regression approach with EA using DEAP (Fortin et al., 2012) , a popular package for symbolic regression research (Quade et al., 2016; Claveria et al., 2016) . We define a set of required primitives: +, -, * , /, x and 1. An expression is represented as a tree where all the primitives are nodes. We start with a population of 10 individual trees with the maximum tree height set to 50. The probability of mating two individuals is 0.1, and of mutating an individual is 0.5 (chosen based on a hyperparameter search, see Appendix I). The limit of number of evaluations for a new offspring is set to 500 so that it is comparable to the number of simulations in MCTS.

H GRAMMAR VARIATIONAL AUTOENCODER:

We would like to point out that the dataset in this paper is more challenging for GVAE than the dataset used in Grammar Variational Autoencoder (GVAE) paper (Kusner et al., 2017) , although it is constructed with fewer production rules. First, the maximum length of production rule sequences in GVAE paper is 15, while the maximum length is 31 in our training set and 43 in holdout set. The RNN decoder usually has difficulties in learning longer sequence due to e.g. exposure bias (Bengio et al., 2015) . Second, while our maximum length is more than doubled, our training set only contains 28837 expressions compared to 100000 in GVAE paper. Samples in our training set are sparser in the syntactic space than those of the GVAE paper.

We trained a GVAE using the the open-sourced code2 with the following modifications: 1) using context-free grammar in our paper 2) setting the max length to 100 so it is comparable with MCTS and EA experiments in our paper. The GVAE paper reported reconstruction accuracy 0.537. On our validation set, the reconstruction accuracy is 0.102.

During each iteration in Bayesian optimization, 500 data in the training set are randomly selected to model the posterior distribution and the model will suggest 50 new expressions. These 50 new expressions will be evaluated and concatenated into the training set. For each symbolic regression task, we ran 5 iterations. The only difference from the setting in Kusner et al. (2017) is that instead of averaging over 10 repetitions, we average over 2 repetitions so the total number of evaluation is 2 × 5 × 50 = 500 for searching each expression, which is comparable to the setting of MCTS and EA.

Due to the restriction on the computational time of Bayesian optimization for GVAE, we evaluate GVAE on a subset of holdout sets. 30 expressions are randomly selected from each of the holdout We also select expressions solved by EA + PW but unsolved by NG-MCTS in Table 1 as examples of symbolic regression results. -200, 200] for the four subplots in the second row to show the discrepancy of expressions on two different scales.

K SYMBOLIC REGRESSION WITH NOISE:

In the main text, the training points D train from the desired expression is noise-free. However, in realistic applications of symbolic regression, measurement noise usually exists. We add a random Gaussian noise with standard deviation 0.5 to expression evaluations on D train and compute ∆g int. and ∆g ext. using evaluations without noise. The results are summarized in Table K .1. The performance of all the methods is worse than that in the noise-free experiments, but the relative relationship 

L SYNTACTIC NOVELTY EXAMPLES:

For ease of presentation, we show syntactic novelties generated by NN that only have one semantically identical expression (i.e., the expression that shares the same simplified expression) in the training set. By comparing each syntactic novelty and its semantically identical expression in the training set (shown in Table L .1), we can observe that the model generates some nontrivial syntactic novelties. For -200, 200] for the four subplots in the second row to show the discrepancy of expressions on two different scales. -200, 200] for the four subplots in the second row to show the discrepancy of expressions on two different scales. ), off (×) and not available (-). If the RMSEs of the best found expression g(x) in interpolation and extrapolation are both smaller than 10 -9 and ∆P [g] = 0, it is solved. If g(x) is non-terminal or ∞, it is invalid. Hard includes expressions in the holdout set which are not solved by any of the six methods. The medians of ∆g train , ∆g int. , ∆g ext. and the median absolute errors of leading powers ∆P [g] for hard expressions are reported. 

:

Model implemented in TensorFlow(Abadi et al.,

2016) and available in submitted materials.

https://github.com/mkusner/grammarVAE

Figures and Tables:

Figure 1: Overview. (a) Exemplary expression. (b) Leading powers of 1/(x + 1) at 0 and ∞. (c) Parse tree of 1/(x + 1). (d) Production rule sequence, the preorder traversal of production rules in the parse tree. (e) Architecture of the model to predict the next production rule from the partial sequence conditioned on desired leading powers. (f) Using (e) to guide MCTS.

The training expressions are sampled sparsely, which are only 10 -23 of the expressions within 31 production rules. Symbolic regression tasks are evaluated on 4250 expressions in holdout sets with M [f ] ≤ 4, = 5, = 6. The challenges are: 1) conditions M [f ] = 5, 6 do not exist in training set; 2) the search spaces of M [f ] = 5, 6 are 10 7 and 10 11 times larger than the training set.

Figure 2: Extrapolation errors of symbolic regression methods in holdout set M [f ] = 5. Each expression is a point, where log 10 ∆g ext. obtained by NG-MCTS is on x-axis and those obtained by other three methods are on the y-axis.

Figure3: Plot of force field expressions found by each method. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . We first discuss the results for holdout set M [f ] ≤ 4. Conventional symbolic regression only fits on the data points D train . EA solves 12.83% expressions, while MCTS only solves 0.54% expressions. This suggests that compared to EA, MCTS is not efficient in searching a large space with limited number of simulations. The median errors ∆P [g] are both 3 for hard expressions (expressions unsolved by all methods), which are large as maximum M [f ] for this set is 4.

Figure 4: Visualizing metrics for conditional production rule generating NN and LHNC (8) on each condition within |P x→0 [f ]| ≤ 9 and |P x→∞ [f ]| ≤ 9. Conditions with M [f ] ≤ 4 are inside the red boundary and points with 0 value are left blank.

each pair of integer leading powers satisfying M [f ] ≤ 4, 1000 shortest expressions are selected to obtain 41000 expressions in total. They are randomly split into three sets. The first two are training (28837) and validation (4095) sets. For the remaining expressions, 50 expressions with unique simplified expressions are sampled from each condition for M [f ] ≤ 4, to form a holdout set with 2050 expressions. In the same way, we also create a holdout set of 1000 expressions for M [f ] = 5 and 1200 expressions for M [f ] = 6. These conditions do not exist in training and validation sets.

Figure A.1 shows the histogram of lengths of the production rule sequence for expressions in the training set and holdout sets. TableA.1 shows the minimum, median and maximum of the lengths of the production rule sequence. The last column, space size, shows the number of possible expressions within the maximum length of the production rule sequence (including those nonterminal expressions). It is computed recursively as follows. Let N * ,i denote the number of possible expressions with length ≤ i and whose production rule sequences start with symbol * , where * can be O, S and T . Then we have

Figure A.1: Histograms of lengths of the production rule sequence.

Figure D.1: Cumulative counts of unique expressions and unique simplified expressions that satisfy the desired condition among expressions generated by the NN model on various desired conditions. Left: desired condition (c (0) , c (∞) ) = (0, 0). Middle: desired condition (c (0) , c (∞) ) = (5, 5). Right: desired condition (c (0) , c (∞) ) = (-2, -3).

Figure E.1: Metrics for NN, NNNC and Random models on each condition within |P x→0 [f ]| ≤ 9 and |P x→inf [f ]| ≤ 9. Each column corresponds to a metric: success rate, mean L1-distance, syntactic novelty rate and semantic novelty rate, respectively. Conditions with M [f ] ≤ 4 are inside the red boundary. Conditions with zero success rate, syntactic novelty rate or semantic novelty rate are left blank in the corresponding plots.

Figure E.2: Metrics for FH and FHNC models on each condition within |P x→0 [f ]| ≤ 9 and |P x→inf [f ]| ≤ 9. Each column corresponds to a metric: success rate, mean L1-distance, syntactic novelty rate and semantic novelty rate, respectively. Conditions with M [f ] ≤ 4 are inside the red boundary. Conditions with zero success rate, syntactic novelty rate or semantic novelty rate are left blank in the corresponding plots.

Figure E.3: Metrics for LH models with different history length on each condition within |P x→0 [f ]| ≤ 9 and |P x→inf [f ]| ≤ 9. Each column corresponds to a metric: success rate, mean L1-distance, syntactic novelty rate and semantic novelty rate, respectively. Conditions with M [f ] ≤ 4 are inside the red boundary. Conditions with zero success rate, syntactic novelty rate or semantic novelty rate are left blank in the corresponding plots.

Figure E.4: Metrics for LHNC models with different history length on each condition within |P x→0 [f ]| ≤ 9 and |P x→inf [f ]| ≤ 9. Each column corresponds to a metric: success rate, mean L1-distance, syntactic novelty rate and semantic novelty rate, respectively. Conditions with M [f ] ≤ 4 are inside the red boundary. Conditions with zero success rate, syntactic novelty rate or semantic novelty rate are left blank in the corresponding plots.

Figure J.1, Figure J.2 and Figure J.3 show eight expressions with M [f ] ≤ 4, M [f ] = 5 and M [f ] = 6, respectively. The symbolic expressions, leading powers, interpolation errors and extrapolation errors of these 24 desired expressions f (x), as well as their corresponding best expressions found by NG-MCTS, denoted by g NG-MCTS (x), and by EA + PW, denoted by g EA+PW (x), are listed in Table J.1.

Figure J.4, Figure J.5 and Figure J.6 show eight expressions with M [f ] ≤ 4, M [f ] = 5 and M [f ] = 6, respectively. The symbolic expressions, leading powers, interpolation errors and extrapolation errors of these 24 desired expressions f (x), as well as their corresponding best expressions found by NG-MCTS and EA + PW are listed in Table J.2.

Figure J.1: Examples of expressions solved by NG-MCTS but unsolved by EA + PW with M [f ] ≤ 4. Each subplot of (a)-(h) demonstrates an expression solved by NG-MCTS but unsolved by EA + PW. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Figure J.2: Examples of expressions solved by NG-MCTS but unsolved by EA + PW with M [f ] = 5. Each subplot of (a)-(h) demonstrates an expression solved by NG-MCTS but unsolved by EA + PW. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Figure J.3: Examples of expressions solved by NG-MCTS but unsolved by EA + PW with M [f ] = 6. Each subplot of (a)-(h) demonstrates an expression solved by NG-MCTS but unsolved by EA + PW. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Figure J.4: Examples of expressions solved by EA + PW but unsolved by NG-MCTS with M [f ] ≤ 4. Each subplot of (a)-(h) demonstrates an expression solved by EA + PW but unsolved by NG-MCTS. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Figure J.5: Examples of expressions solved by EA + PW but unsolved by NG-MCTS with M [f ] = 5. Each subplot of (a)-(h) demonstrates an expression solved by EA + PW but unsolved by NG-MCTS. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Figure J.6: Examples of expressions solved by EA + PW but unsolved by NG-MCTS with M [f ] = 6. Each subplot of (a)-(h) demonstrates an expression solved by EA + PW but unsolved by NG-MCTS. Grey area is the region to compute interpolation error ∆g int. and light blue area is the region to compute extrapolation error ∆g ext. . The display range of y-axis is [-5, 5] for the four subplots in the first row and [-200, 200]  for the four subplots in the second row to show the discrepancy of expressions on two different scales.

Results of symbolic regression methods with noise. Search expressions in holdout sets M [f ] ≤ 4, M [f ] = 5 and M [f ] = 6 with data points on D train and / or leading powers P x→0 [f ] and P x→∞ [f ]. The options are marked by on ( √

x)/((1 + x) + x) A/(B + B) = A/B/2 1 -x/((1 + x) + (1 + x)) 1 -(x/(1 + x))/(1 + 1) A + B = B + A x/((x/(1 + x)) + x + x) x/(x + (x/(1 + x)) + x) A/(B * C) = A/B/C (1 -(x/(1 -x)))/(x + x * x) (1 -(x/(1 -x)))/(x + 1)/x A/C + B/C = (A + B)/C (1 + (1/x))/(x -(x * x) -1) ((1 + x)/x)/(x -1 -(x * x))



Table content: <table><tr><td>{•} | on training points D train :</td></tr><tr><td>{1.2, 1.6, 2.0, 2.4, 2.8}, points in interpolation region D interpolation : {1.4, 1.8, 2.2, 2.6} and points</td></tr><tr><td>in extrapolation region D extrapolation : {5, 6, 7, 8, 9}.</td></tr></table>

Results of symbolic regression methods. Search expressions in holdout sets M [f ] ≤ 4, M [f ] = 5 and M [f ] = 6 with data points on D train and / or leading powers P x→0 [f ] and P x→∞ [f ].The options are marked by on ( √ ), off (×) and not available (-). If the RMSEs of the best found expression g(x) in interpolation and extrapolation are both smaller than 10 -9 and ∆P [g] = 0, it is solved. If g(x) is non-terminal or ∞, it is invalid. Hard includes expressions in the holdout set which are not solved by any of the six methods. The medians of ∆g train , ∆g int. , ∆g ext. and the median absolute errors of leading powers ∆P [g] for hard expressions are reported. Evaluated on the subset of holdout sets. Hard are unsolved expressions in the subset (Appendix H).

Table content: <table><tr><td colspan="2">M [f ] Method</td><td>Neural</td><td colspan="2">Objective Function</td><td>Solved</td><td>Invalid</td><td/><td/><td>Hard</td><td/><td/></tr><tr><td/><td/><td colspan="5">Guided Dtrain Px→0,∞[f ] Percent Percent</td><td>Percent</td><td colspan="4">∆gtrain ∆gint. ∆gext. ∆P [g]</td></tr><tr><td>≤ 4</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td colspan="2">0.54% 0.24% 0.20% 71.22% 0.00% 2.93% 0.00% 0.39% 12.83% 3.32% 23.37% 0.44%</td><td>23.66%</td><td>0.728 -0.967 0.225 0.186 0.376</td><td>0.598 2.069 0.836 0.194 0.162 0.322</td><td>0.723 2.823 0.541 0.084 0.358 0.152</td><td>3 1 2 0 3 0</td></tr><tr><td/><td>GVAE  ‡ GVAE + PW  ‡</td><td>--</td><td>√ √</td><td>× √</td><td>10.00% 10.00%</td><td>0.00% 0.00%</td><td>90.00% 90.00%</td><td>0.217 0.386</td><td>0.159 0.324</td><td>0.599 0.056</td><td>2 0</td></tr><tr><td>= 5</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td colspan="2">0.00% 0.00% 0.10% 32.10% 0.00% 3.90% 0.00% 3.50% 2.90% 4.20% 9.20% 2.30%</td><td>58.00%</td><td>0.857 -1.105 0.247 0.227 0.366</td><td>0.738 1.890 0.914 0.229 0.204 0.365</td><td>0.950 1.027 0.600 0.020 0.155 0.109</td><td>5 3 4 0 4 0</td></tr><tr><td/><td>GVAE  ‡ GVAE + PW  ‡</td><td>--</td><td>√ √</td><td>× √</td><td>0.00% 3.33%</td><td>0.00% 0.00%</td><td>100.00% 96.67%</td><td>0.233 0.649</td><td>0.259 0.565</td><td>0.164 0.383</td><td>4 2</td></tr><tr><td>= 6</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td>0.00% 0.00% 0.08% 17.33% 1.25% 4.92%</td><td>6.33% 0.00% 7.33% 0.17% 5.08% 6.92%</td><td>71.25%</td><td>1.027 -1.228 0.236 0.219 0.329</td><td>0.819 2.223 1.051 0.209 0.191 0.285</td><td>0.852 7.145 0.891 0.008 0.084 0.047</td><td>6 4 4 0 5 0</td></tr><tr><td/><td>GVAE  ‡ GVAE + PW  ‡</td><td>--</td><td>√ √</td><td>× √</td><td>0.00% 0.00%</td><td>0.00% 0.00%</td><td>100.00% 100.00%</td><td>0.260 0.595</td><td>0.206 0.436</td><td>0.037 0.087</td><td>5 3</td></tr><tr><td>‡</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

The upper right cluster in each plot represents expressions unsolved by both methods. Most of the plotted points are above the 1:1 line (dashed), which shows that NG-MCTS outperforms the others for most unsolved expressions in extrapolation. Examples of expressions solved by NG-MCTS but unsolved by EA + PW and vice versa are presented in Appendix J. We also perform similar experiments with Gaussian noise on D train and NG-MCTS still outperforms all other methods ( Appendix K). Results of force field expressions found by each method.

Table content: <table><tr><td>Method</td><td>Expression Found</td><td>∆gtrain ∆gint. ∆gext. ∆P [g]</td></tr></table>

Metrics for conditional production rule generating NN and baseline models.

Table content: <table/>

FH behaves as a look-up table based method (i.e., sampling from the training set) so it has 100% success rate in in-sample condition region. However, it is not able to generate any novel expressions. NN has the best performance on the syntactic and semantic novelty rates in both the insample (35% and 2.7%) and out-of-sample (1416 and 1084) condition regions by a significant margin. This indicates the generalization ability of the model to generate unseen expressions matching a desired condition. It is worth pointing out that NNNC performs much worse than NN, which indicates that NN is not only learning the distribution of expressions in the dataset, but instead is also learning a conditional distribution to map leading powers to the corresponding expression distributions.

Table content: <table><tr><td/><td/><td>Sem (%)</td><td colspan="4">Total Num Expressions Num Conditions Suc Syn Sem with Suc</td><td/><td colspan="2">Mean L1-Dist</td><td/></tr><tr><td>NN</td><td>35.0</td><td>2.7</td><td colspan="2">1465 1416</td><td>1084</td><td>115</td><td>0.8</td><td>1.4</td><td>2.5</td><td>4.3</td></tr><tr><td>NNNC</td><td>1.8</td><td>0.2</td><td>7</td><td>7</td><td>7</td><td>7</td><td>4.0</td><td>5.6</td><td>6.5</td><td>7.5</td></tr><tr><td>Random</td><td>0.7</td><td>0.0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>10.9</td><td colspan="3">11.7 12.4 12.6</td></tr><tr><td>FH</td><td>0.0</td><td>0.0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0.0</td><td colspan="3">18.0 18.0 18.0</td></tr><tr><td>FHNC</td><td>0.0</td><td>0.0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>4.2</td><td>5.7</td><td>6.6</td><td>7.5</td></tr><tr><td>LH (2)</td><td>5.0</td><td>1.1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>3.1</td><td colspan="3">18.0 18.0 18.0</td></tr><tr><td>LH (4)</td><td>10.3</td><td>2.1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>2.5</td><td colspan="3">18.0 18.0 18.0</td></tr><tr><td>LH (8)</td><td>14.6</td><td>2.4</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1.8</td><td colspan="3">18.0 18.0 18.0</td></tr><tr><td>LH (16)</td><td>1.2</td><td>0.1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0.1</td><td colspan="3">18.0 18.0 18.0</td></tr><tr><td>LHNC (2)</td><td>1.4</td><td>0.3</td><td>7</td><td>7</td><td>7</td><td>6</td><td>4.2</td><td>5.6</td><td>6.9</td><td>7.5</td></tr><tr><td>LHNC (4)</td><td>1.2</td><td>0.3</td><td>6</td><td>6</td><td>6</td><td>6</td><td>3.9</td><td>5.7</td><td>6.4</td><td>7.3</td></tr><tr><td>LHNC (8)</td><td>1.5</td><td>0.3</td><td>8</td><td>8</td><td>8</td><td>6</td><td>4.2</td><td>5.9</td><td>6.7</td><td>7.6</td></tr><tr><td>LHNC (16)</td><td>0.2</td><td>0.1</td><td>3</td><td>3</td><td>3</td><td>2</td><td>4.3</td><td>5.7</td><td>6.6</td><td>7.5</td></tr><tr><td>example,</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Examples of top-3 probable expression completions for different desired conditions.

Table content: <table/>

1: Minimum, median and maximum of the lengths of the production rule sequence. Space size is the number of possible expressions within the maximum length of the production rule sequence.

Table content: <table><tr><td>Name</td><td/><td>Length</td><td/><td colspan="2">Number of Expressions Space Size</td></tr><tr><td/><td colspan="3">Min Median Max</td><td/><td/></tr><tr><td>training M [f ] ≤ 4</td><td>3</td><td>19</td><td>31</td><td>28837</td><td>2.2 × 10 27</td></tr><tr><td>holdout M [f ] ≤ 4</td><td>7</td><td>19</td><td>31</td><td>2050</td><td>2.2 × 10 27</td></tr><tr><td>holdout M [f ] = 5</td><td>15</td><td>27</td><td>39</td><td>1000</td><td>8.9 × 10 34</td></tr><tr><td>holdout M [f ] = 6</td><td>11</td><td>31</td><td>43</td><td>1200</td><td>5.8 × 10 38</td></tr></table>

1: Comparison with GVAE on subset of holdout sets. this section, we select expressions solved by NG-MCTS but unsolved by EA + PW in Table1as examples of symbolic regression results.

Table content: <table><tr><td colspan="2">M [f ] Method</td><td>Neural</td><td colspan="2">Objective Function</td><td>Solved</td><td>Invalid</td><td/><td/><td>Unsolved</td><td/></tr><tr><td/><td/><td colspan="5">Guided Dtrain Px→0,∞[f ] Percent Percent</td><td>Percent</td><td colspan="4">∆gtrain ∆gint. ∆gext. ∆P [g]</td></tr><tr><td>≤ 4</td><td>NG-MCTS GVAE GVAE + PW</td><td>√ --</td><td>√ √ √</td><td>√ × √</td><td colspan="2">83.33% 0.00% 10.00% 0.00% 10.00% 0.00%</td><td>16.67% 90.00% 90.00%</td><td>0.436 0.217 0.386</td><td>1.000 0.159 0.324</td><td>0.009 0.599 0.056</td><td>0 2 0</td></tr><tr><td>= 5</td><td>NG-MCTS GVAE GVAE + PW</td><td>√ --</td><td>√ √ √</td><td>√ × √</td><td colspan="2">60.00% 0.00% 0.00% 0.00% 3.33% 0.00%</td><td>40.00% 100.00% 96.67%</td><td>0.088 0.233 0.649</td><td>0.086 0.259 0.565</td><td>0.011 0.164 0.383</td><td>0 4 2</td></tr><tr><td>= 6</td><td>NG-MCTS GVAE GVAE + PW</td><td>√ --</td><td>√ √ √</td><td>√ × √</td><td colspan="2">10.00% 0.00% 0.00% 0.00% 0.00% 0.00%</td><td>90.00% 100.00% 100.00%</td><td>0.306 0.260 0.595</td><td>0.266 0.206 0.436</td><td>0.009 0.037 0.087</td><td>0 5 3</td></tr><tr><td colspan="8">J EXAMPLES OF SYMBOLIC REGRESSION RESULTS</td><td/><td/><td/></tr></table>

2: Examples of expressions solved by EA + PW but unsolved by NG-MCTS. This table shows the desired expressions f (x) and their corresponding best expressions found by NG-MCTS g NG-MCTS (x) and by EA + PW g EA+PW (x). The leading powers P x→0 [•], P x→∞ [•], interpolation error ∆g int. and extrapolation error ∆g ext. are reported for each expression.

Table content: <table><tr><td>∆gext.</td></tr><tr><td>Px→0[•] Px→∞[•] ∆gint.</td></tr><tr><td>Expression</td></tr></table>

Table L.1: Examples of syntactic novelties to demonstrate what the NN model learned. Semidentical expression refers to the expression in the training set that shares the same simplified expression as the corresponding syntactic novelty. The first column shows the mathematical rules a human needs to know and apply to derive each syntactic novelty from its semantically identical expression in the training set.

Table content: <table><tr><td colspan="2">M [f ] Method</td><td>Neural</td><td colspan="2">Objective Function</td><td>Solved</td><td>Invalid</td><td/><td/><td>Hard</td><td/><td/></tr><tr><td/><td/><td colspan="10">Guided Dtrain Px→0,∞[f ] Percent Percent Percent ∆gtrain ∆gint. ∆gext. ∆P [g]</td></tr><tr><td>≤ 4</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td>0.39% 0.34% 0.34% 24.29% 4.44% 4.34%</td><td>0.54% 0.00% 0.49% 0.05% 3.95% 0.39%</td><td>71.41%</td><td>0.689 -0.887 0.432 0.399 0.385</td><td>0.351 1.003 0.643 0.449 0.223 0.489</td><td>0.371 0.865 0.825 0.256 0.591 0.260</td><td>3 1 2 0 3 0</td></tr><tr><td>= 5</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td>0.00% 0.00% 0.00% 10.80% 1.00% 1.80%</td><td>4.00% 0.00% 3.00% 0.10% 3.10% 1.80%</td><td>83.40%</td><td>0.931 -0.944 0.558 0.480 0.448</td><td>0.599 1.394 0.817 0.430 0.256 0.382</td><td>0.972 1.000 0.727 0.103 0.266 0.122</td><td>5 3 5 0 4 0</td></tr><tr><td>= 6</td><td>MCTS MCTS (PW-only) MCTS + PW NG-MCTS EA EA + PW</td><td>× × × √ --</td><td>√ × √ √ √ √</td><td>× √ √ √ × √</td><td>0.00% 0.00% 0.00% 10.33% 0.67% 1.42%</td><td>9.75% 0.00% 8.83% 0.08% 4.58% 6.50%</td><td>78.75%</td><td>0.960 -1.122 0.426 0.463 0.388</td><td>0.762 1.024 0.807 0.205 0.427 0.369</td><td>0.888 0.861 0.163 0.009 0.852 0.065</td><td>6 4 4 0 5 0</td></tr></table>

Peer review (iclr2020):
Summary: The authors of this paper propose a novel approach for symbolic regression. The simulation results demonstrate that the proposed approach can find a better function g producing similar results as desired function f by providing the additional information – the leading power of function f. Paper strength: 1. The proposed NG-MCTS is elegant to solve the problem of symbolic regression. 2. Experiment results illustrate the superiority of the proposed approach Paper weakness: 1. I can follow most of the mathematics in the paper. But the most confusing part for me is why you feed the random partial sequence for training. Besides, how you do inference to generate a sequence of the production rules. 2. What is the final objective function? If I do not misunderstand, it could be the cross-entropy loss between the output of GRU and the next target production rule, RMSE and the error on leading power. Then how you optimize it? Please describe more details about this. 3. The authors of this paper introduce more information – leading power of desired function for symbolic regression but they incorporate the additional information by introducing a simple loss function term. How about the performance of baseline approaches with those kinds of information? 4. The whole systems seem like very complicate and it would be more interesting if the authors provide sufficient ablations to decompose the complex algorithm.