{'fraction_non_alphanumeric': 0.03729718339537824, 'fraction_numerical': 0.01643604691999719, 'mean_word_length': 4.704821690930947, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Reinforcement Learning Approach to Domain-Knowledge Inclusion Using Grammar Guided Symbolic Regression

In recent years, symbolic regression has been of wide interest to provide an interpretable symbolic representation of potentially large data relationships. Initially circled to genetic algorithms, symbolic regression methods now include a variety of Deep Learning based alternatives. However, these methods still do not generalize well to real-world data, mainly because they hardly include domain knowledge nor consider physical relationships between variables such as known equations and units. Regarding these issues, we propose a Reinforcement-Based Grammar-Guided Symbolic Regression (RBG2-SR) method that constrains the representational space with domain-knowledge using context-free grammar as reinforcement action space. We detail a Partially-Observable Markov Decision Process (POMDP) modeling of the problem and benchmark our approach against state-of-the-art methods. We also analyze the POMDP state definition and propose a physical equation search use case on which we compare our approach to grammar-based and non-grammarbased symbolic regression methods. The experiment results show that our method is competitive against other state-of-the-art methods on the benchmarks and offers the best error-complexity trade-off, highlighting the interest of using a grammar-based method in a real-world scenario.


Introduction
Finding a generic symbolic representation from an observation set has been of a long interest in the Physics community. Long before Newton's fallen apple, scientists have been working on discovering symbolic relationships that satisfy observations, current knowledge, and already known equations. In astronomy, for example, the search for a closed-form solution has been an ongoing interest notably for Kepler's study of planetary motion and still today for the modeling of albedos (Heng et al, 2021). However, in recent years, the number of observations available has surged thanks to the increase of sensor monitoring systems in Big Data environments and, finding the best possible equation to model complex systems is now a tedious task without computer assistance. This automatic equation search task is known as Symbolic Regression (SR). It aims at finding a symbolic function f that matches the relationship f (X) = y between an observation set X ∈ R n described with n variables and a target variable y ∈ R to explain y from X. SR was, until recently, mainly performed using Genetic Programming (GP) (Koza, 1990), a family of techniques that draws inspiration from Darwinian evolution to search for solutions automatically. However, the initial GP formulation is not suited for large search spaces (Ebner, 1999). To cope with this limitation, since the early works on GP (Koza, 1990), multiple extensions have been proposed to restrict the search space, such as Strongly Typed Genetic Programming (Montana, 1995) or Grammar Guided Genetic Programming (G3P) (Koza et al, 2006). These methods share the interesting property of providing a way to enforce domain knowledge into the learning process (Ratle and Sebag, 2000) and are thus applied to various real-world problems (Cherrier et al, 2019a).
Meanwhile, Deep Learning (DL) methods have progressively become ubiquitous because of their high representational capacity when trained on large datasets (Russakovsky et al, 2015). However, in the general case, DL often has a "black box" behavior and lack of interpretability. Thus, combining DL and SR by taking advantage of the computational capacity of DL and the expressiveness of SR would help provide more human-readable results. Along with the rapid development of DL, there have also been a growing interest in Reinforcement Learning (RL) methods to solve decision-making problems. RL is a Machine Learning paradigm inspired by behavioral psychology concerned with how an agent should act in an environment to maximize a cumulative reward (Sutton and Barto, 2018). Deep reinforcement learning (Deep RL) approaches especially show great power in solving complex sequential problems. Recent work in these both domains also offers viable alternatives to GP on SR tasks Udrescu and Tegmark, 2020). Still, unlike GPbased methods, they are not yet able to include sophisticated domain-related knowledge or constraints.
In this work, we propose a method called Reinforcement Based Grammar Guided Symbolic Regression (RBG2-SR) to tackle SR with a Deep Reinforcement Learning (Deep RL) approach using a Backus-Naur Form (BNF) (Knuth, 1964) Context Free Grammar action space which constrains the solution space only to domain-viable solutions. Our method enforces domain-related constraints as grammatical rules in a human-readable manner. These constraints are then translated into directly interpretable symbolic outputs. We also propose a general Partially-Observable Markov Decision Process (POMDP) (Kaelbling et al, 1998) modeling of the SR problem. We show that our approach performs better than other tested methods with the same grammar on the tested benchmarks. We also offer comparative results on a real-world scenario to show how our approach could use a BNF grammar to take advantage of expert knowledge.
The rest of this paper is organized as follows. First, Section 2 summarizes related state-of-the-art works. In Sections 3 and 4, we respectively describe the proposed method and its corresponding experimental results. Finally, Section 5 offers concluding remarks and perspectives.

Symbolic Regression
Symbolic Regression (SR) is the process of searching for a symbolic relationship, also called expression, that accurately matches a given dataset. Found expressions can be represented as trees where nodes contain operations and variables. SR was early investigated with Genetic Programming (GP) approaches (Koza, 1990(Koza, , 1992. They iteratively evolve a population of individuals (each individual representing an expression) across multiple generations through evolutionary operations. Many works today are descendants of these works, and propose to overcome the problems of GP, such as bloat (Silva, 2008) (the individual's complexity explosion over the generations) or convergence (Rosca, 1996) issues, for example, with multi-objective strategies (Tamaki et al, 1996) or a partial derivative-based error fitness (Schmidt and Lipson, 2009).
Driven by the current need for more interpretable models, non-GP-based methods have been developed to tackle SR. These methods propose to take advantage of the computational capacity of neural networks (Hornik et al, 1989) while providing an interpretable solution. For instance, they offer to encode the expression in the neural network structure and activation functions (Sahoo et al, 2018;Kim et al, 2020), to predict a string expression (Anjum et al, 2019) with Recurrent Neural Networks, or to use Deep Reinforcement Learning (Deep RL) as a search engine . By combining partial derivative and neural networks, AI Feynman method (Udrescu and Tegmark, 2020) propose to make use of simplifying properties (such as units, symmetry, separability.. etc.) intrinsically present in physical expressions to repeatedly cut the global SR problem into simpler ones with fewer variables. Other methods that do not use neural networks rely on bayesian optimization (Jin et al, 2019) or perform nonlinear basis function expansion (McConaghy, 2011).
However, as designed, these methods do not insert custom knowledge and expertise into the symbolic expression construction.

Knowledge Insertion by Constraints
Because the search space in SR is very large and can leads to local optima, it is relevant to restrict the function space by removing sub-optimal search regions. In line with Koza's work on GP, a preliminary solution called Strongly Typed GP (STGP) (Montana, 1995) proposed to enforce data types constraints for computer program search, which decreases the search time and improves the generalizability of the found solutions. Regarding SR for physical laws (re)discovery, other kinds of constraints can be considered, such as physical units. More precisely, as each variable comes with its physical unit, arbitrarily combined variables can produce illegal unit combinations. To this end, dimensionally awareness GP (Keijzer and Babovic, 1999) was initially proposed to take into account unit knowledge in GP by minimizing the distance to a legal unit in the fitness function. Note that, more than just constraining the search space, these constraints enforce structured knowledge and expertise about the problem within the learning. This knowledge can also take the form of ontologies (Prieschl et al, 2019) to include prior knowledge as additional input features. However, both dimensionally-aware GP (Keijzer and Babovic, 1999) and ontology-guided GP (Prieschl et al, 2019) do not ensure to produce dimensionally valid expressions.
Toward this goal, the definition of explicit constraints can guarantee to produce only legal expressions with respect to the constraints. These constraints can be grammatical (Whigham et al, 1995) or ontological (Lucena-Sánchez et al, 2021). In Grammar-Guided Genetic Programming (G3P) (Whigham et al, 1995), also called Grammar-Based GP, a Context-Free Grammar (CFG) (Cremers and Ginsburg, 1975) is used to define constraint rules. Grammatical rules allow defining physical units, thanks to which G3P has found a variety of industrial applications (Crochepierre et al, 2021;Cherrier et al, 2019b). CFGs are often written in Backus-Naur form (BNF) (Knuth, 1964), which is made of: • Non-Terminals also called symbols (ex. <s>). One of the non-terminals is called start symbol. • Terminals, character strings to replace non-terminals (ex. "a", "b", "+").
The set of terminals can contain input features and operators to combine features. • Rules which defines how the terminals and non-terminals are connected.
A set of rules for a specific symbol is called a production rule where rules are separated by a vertical line |, and ::= means "defined as" (for example <s> ::= "a" | "b" | <s>+<s>). Eventually, a grammar contains multiple production rules, one per symbol. To select which rule will replace a given symbol <symbol>, uniform sampling is made in <symbol> production. McKay et al. (McKay et al, 2010) provides a detailed survey of G3P strategies and genetic operations. Because of its ability to restrict the symbol space, CFG structures have been a privileged topic of study. Outside SR, CFGs found other applications, such as in Grammar Variational Autoencoder (Kusner et al, 2017) where they are used for symbolic data representation along with molecule representation.
Nowadays, probabilistic CFG (noted PCFG) (Sakakibara, 2017) are prefered as they allow to weigh the importance of a rule. In addition to CFG rules, they assign a probability to each rule in a production rule so that all probabilities for a given symbol add up to 1. At the end of a production rule, a list preceded by a keyword probs and a double vertical line defines the probabilities associated to the rules. The production rule structure now becomes: <symbol> ::= rule1 | rule2 |... || probs [prob_r1,prob_r2,...] PCFGs are of particular importance as they allow estimating the probabilities associated with each rule and discard unuseful rules. Probability distributions can be updated according to sampled expressions using Linear Genetic Programming (Sotto and de Melo, 2017) or Monte Carlo sampling (Brence et al, 2021).

Reinforcement Learning
Reinforcement Learning (RL) is a Machine Learning approach to solving Markov Decision Processes (MDPs), where the MDP is mainly defined by its state space, action space, state transitions probabilities, and reward. In the RL paradigm (Sutton and Barto, 2018), an agent learns to achieve a task by interacting with its environment at discrete time steps. At each time step, the agent chooses an action among available actions according to a given policy. Next, the action is sent to the environment, which gives back a reward feedback and a new state to the agent. Then, the agent can learn from the received reward signals to improve its policy. RL strategies are mainly valued-based like Deep Q-Network (Mnih et al, 2015), or policy-based like REINFORCE algorithm (Williams, 1992). Value-based RL learns a value function and deduce a policy from values. In contrast, policy-based RL explicitly learns a policy π and keeps it in memory during learning. In Deep Reinforcement Learning (Deep RL), value and policy functions are approximated using neural networks. For example, the REINFORCE algorithm (Williams, 1992) at its core uses the policy gradient theorem to update the probability distribution of actions. Actor-Critic algorithms (Konda and Tsitsiklis, 2000) are inheriting from both strategies by trying to learn alongside a policy and its value to reduce variance and improve converge.
RL has been applied to a variety of domains in pattern recognition (Piñol et al, 2012;Khurana et al, 2018;Bertsekas, 2019), from feature construction (Khurana et al, 2018) to feature-based aggregation (Bertsekas, 2019). Including knowledge in RL is of paticular importance especially for safety issues (Alshiekh et al, 2018) where shielding strategy is used to correct actions if the chosen one causes a violation of some specified sort. However, few works have focused on building a symbolic representation for data. Recent work on this topic includes the use of RL to search among a library of operators and features for SR  or the creation of symbolic computer programs (Verma et al, 2018).
Regarding interpretability, Deep Learning and RL approaches mostly produce complex solutions that are often hard to interpret. It is especially true in environments where most agents perform according to a black-box policy. To make these black-box approaches more grey and learn more interpretable RL policies, recent work proposes to learn symbolic policies either with GP (Hein et al, 2018) or Deep RL . However, even if they have interpretable outputs, most of these solutions do not include domain-related knowledge within the learning to ensure that the output solutions follow the prior knowledge.

Definition of the Reinforcement Learning Environment
In this work, we adopt a Markov Decision Process (MDP) (Sutton and Barto, 2018) modeling of the SR problem. More precisely, we choose to consider a Partially-Observable Markov Decision Process (POMDP) (Kaelbling et al, 1998) in a finite episodic setting with maximum horizon H. This section is dedicated to the definition of the main components of the MDP in a RL setting, namely: state and action spaces and reward. In Section 3.1.1, we define State and action spaces for SR and Section 3.1.2 details the reward definition along with its properties. Section 3.1.3 combine the definitions from Sections 3.1.1 and 3.1.2 to define the whole POMDP. Finally, Section 3.1.4 details how to learn a policy over action, used to generate the action probabilities given the current state.

Grammatical state and action space
Let us consider a RL setting where the overall task is to find an optimal symbolic function f * so that f * = argmin f ∈F G ||y − f (X)||, with F G being the function space accessible from a given grammar G. The grammar is defined by the tuple (σ start , (σ nt ) nt∈N T , (σ t ) t∈T , ρ, Ψ) with σ start the start symbol of the grammar, (σ nt ) nt∈N T a set of non-terminals, (σ t ) t∈T a set of terminals, ρ the production rules to combine terminals/non terminals, and Ψ the probability associated to each production rule. Figures 1 and 4 in Section 4 show example BNF Grammar inspired by the work of Sotto and Melo (Sotto and de Melo, 2017). Given this notation, we propose to define the construction of f * as a sequential decision making problem where an agent sequentially chooses rules in the grammar to build up the function f . In this grammatical space, we first specify the maximal number of steps to create f , called the maximal horizon H. We then define for each step h ∈ [0, ..., H] an action a h as the selection of a rule in the production rule accessible from the current state. We also propose to define the state s h at step h by .., a h−1 ] all previously selected action, a parent h the action taken by the parent in the parse tree, a siblings h the action taken by each already computed siblings in the parse tree, d h the depth of the expression tree at step h, σ h the current type of symbol to find at step h, m h a mask over accessible actions from the current state symbol σ h and η h hidden information about the current state. We call sibling nodes those nodes with the same depth as the current node, and parent node the node above the current node in the parse tree. Alternative state definitions will be tested in Section 4.1.3. Initially, we define the state by the start symbol σ start , previously selected actions is an empty list, m 0 masks out inaccessible actions from the initial symbol, and hidden information η 0 is randomly initialized. A trajectory of actions is associated to each symbolic function f k . We consider a case where the grammar is chosen and constructed so that there is always at least one action accessible at each step.
Algorithm 1 Single episode sampling, returns one function f per episode.
queue, a past , a parents , a siblings , f ← empty 3: for h in H do 8: state ← ((a past , a parent , a siblings , d, σ, m, η) 9: action probs, η ← π θ (state) π θ is described in Figure 3(a) 10: action ← sample(action probs) Corresponds to "Action Sampling" blue box in Figure 3 11: a past ← append(a past , action) 12: end for return f 20: end function The Algorithm 1 details the episode sampling procedure. A visual example of this expression sampling is provided in Figure 1. From a given grammar with Fig. 1 "x 9 + x 1 " Expression and trajectory generation from a given grammar. A simplistic grammar is given (left), with actions numbered from 1 to 16 and start symbol < exp >. On the right, a trajectory is sampled from this grammar, from which we define the corresponding parse tree and symbolic expression 16 actions and <exp> as start symbol, we generate a trajectory of actions. The trajectory construction goes as follow: 1. We begin by selecting the first action among accessible actions from the start symbol < exp >: either actions 1 or 2. Given the probabilities associated with these actions, we perform a weighted sampling with the probabilities listed in probs as weights. Action 2 is sampled with the rule "< b >". As this rule contains the non-terminal symbol < b > we need to replace it with a rule from the grammar "accessible" for the symbol < b >. 2. The 3 rd row in the grammar defines actions accessible from < b >: actions 5 or 6. Given the weights, we sample action 5, "< i > + < i >". It contains two non-terminal symbols (< i > and < i >) that need to be replaced in the next steps. 3. The non-terminal symbol replacement is performed on one symbol at a time in a depth-first search manner, by looping over the first symbol until reaching a terminal symbol. We iterate this procedure until either the maximal trajectory length is reached or all non-terminal symbols encountered in all action selections have been replaced by a terminal value.

Reward definition
The standard metric to minimize in SR is the Mean Squared Error (MSE).
To match the RL definitions where the reward function is a monotonically increasing function, we use the squashing function 1 1+x .
As the function f can only be evaluated at the end of the episode when the function is complete, the reward r equals 0 until the final step is reached and values 1 1+M SE(y,ŷ) at step H, whereŷ = f (X) as shown in Equation 1. We also note the expected cumulative reward R and highlight that R = h r h = r H . The sparsity property is used in Section 3.2.3 to simplify the REINFORCE algorithm (Williams, 1992) loss function.

Partially-Observable Markov Decision Process
Given the previous space, action and reward definitions, the POMDP we consider here is defined by a tuple (S, A, r, P, H, Ω, O) with S the state space, A the action space, r : S × A − → [0, 1] the reward function, P : S × A − → [0, 1] the transition kernel, Ω = (o 1 , o 2 , ..., o K ) a set of observations and O a set of conditional observation probabilities O(o|s , a). We write as P (s |s, a) the probability of having a transition to state s ∈ S when taking action a in state s. A POMDP setting is here considered to mitigate the fact that action effects are uncertain until the final state is reached and that the state might be partially observable. Each episode k builds up a trajectory of actions τ k = (a h k ) h∈[H] that ends either when the maximum horizon is reached or when the function f k constructed by τ k is complete and can be evaluated.

Policy optimization
Given the grammatical action space defined in Section 3.1.1, there exists an optimal policy π * capable of generating a trajectory as close as possible to the symbolic function f * . Depending on the choice of a judiciously constructed grammar, it is even possible to generate an optimal trajectory exactly corresponding to f * .
To search for this optimal symbolic function, we propose to learn a policy π θ , parametrized by a vector θ to generate the weights of grammatical action rules at each step h ∈ [H] of the trajectory, from which we sample the next action. More precisely, to build f , the stochastic policy π θ assigns a probability vector to accessible actions from a given state. At each step h, the action is then sampled according to the probability vector given by π θ (s h , o h ). Using the reward defined in Section 3.1.2, we iteratively update the parameters of π θ to sample, in the future, more relevant trajectories with respect to the defined reward. The red elements define the results of the sampling on the grammar. Grammar weights are updated by searching for the target expression "x 9 + x 1 " Using the same grammar example from Figure 1, we show in Figure 2 how this grammar could be updated when trained on the search of the expression "x 9 + x 1 ". After training, all unnecessary grammatical rules now have a low or zero probability, and it is easier to generate the correct target expression with this grammar. Fig. 3 Neural Network architecture to learn π θ (Better seen in color). (a) shows the recurrent architecture to predict an action from a given state as described in Algorithm 1 lines 9 and 10. (b) is the color legend we use in (a) and (c). (c) represents the POMDP agent-environment interaction 3.2 Learning π θ and finding f * : Implementation Details

POMDP Modeling with a Recurrent Neural Network
In order to find an optimal function f * , we propose to use π θ as an exploration tool, trained to emphasize exploration on the most relevant regions of the grammatical space. To do so, we learn π θ with the policy-gradient (PG) algorithm called REINFORCE (Williams, 1992) on the neural network architecture described in Figure 3. In Figure 3(a), we describe (using the legend in Figure  3(b)) the recurrent architecture to predict an action a h from a given state s h at step h. To handle recurrency we choose Long Short-Term Memory (LSTM) (Hochreiter and Schmidhuber, 1997) networks. They are a purposely designed structure to capture long term temporal dependencies and they are also able to model the partial observability of the above-mentioned MDP (Wierstra et al, 2007).
The neural network takes as input the state s h (in orange), and hidden state η h (in grey). All state-input features are encoded either using convolutions or feed-forward layers depending on their shape. The encoded state segments are then concatenated and handed to a LSTM cell along with the hidden state η h . This recurrent cell both outputs an estimation of the observations for the next hidden state η h+1 , and an encoding of all actions in the grammar. This actions encoding is then masked using the state mask m h (see Section 3.2.2) and outputs the distribution over actions accessible from the current state. Then, we sample action s h according to the distribution over accessible action.

Invalid action Masking
To handle grammatical constraints on the action space, we propose to use at each step h a mask over inaccessible actions m h in the current state s h . Similarly to what is done by Huang and Ontañón (2020), invalid actions are masked through element-wise multiplication with a large negative value and followed by softmax function.

Exploration Driven Cost Function
REINFORCE algorithm (Williams, 1992) objective function, simplified with Equation 1, is J θ = E π [R h log π θ (a h |s h )] from which we specify the optimal policy π * θ = argmax θ J θ . As PG methods such as REINFORCE tends to have a large variance (Chung et al, 2021), it is common practice to subtract a baseline to the batch, such as a moving average of rewards across batches. However, SR is only interested in building a policy that maximizes the bestperforming trajectories found during training, and these strategies might have a slow convergence because many sampled trajectories are irrelevant and have a low or zero final reward. Based on these comments, Petersen et al (2021) proposed a risk-seeking policy gradient, which only compute the cost function based on the top− quantile of the expected rewards R , i.e., the most relevant trajectories of the batch: They also added a entropy term H, weighted by λ H , to encourage exploration: Eventually, our final cost function becomes: Given this cost function, the Neural network architecture is trained using the Algorithm 2.

Brief summary of RBG2-SR method
To summarize, our proposed RBG2-SR method performs constrained SR where a grammatical structure restricts symbolic expressions creation. We adopt a POMDP setting with a finite horizon H, for which the associated reward r h is zero until the full expression f has been generated. A grammar defines a set of constraint rules used for the expression construction, which masks the nonaccessible actions at each step h ∈ [H]. The weights of the actions accessible at the current time step are generated by π θ , a neural network learned with the REINFORCE algorithm, using a risk-seeking with entropy term loss.

Experiments
In this section, we describe two experiments. The first one compares the proposed method with several state-of-the-art (Whigham et al, 1995;Sotto and de Melo, 2017) SR solutions on reference benchmarks (Uy et al, 2011;Keijzer, 2003;Vladislavleva et al, 2009;Pagie and Hogeweg, 1997). In the second experiment, we present a feature exploration use case where we show an application of our approach on a real-world dataset with an unknown relationship to uncover. Both data and code for this benchmark are freely available on Github 1 .

Benchmarked methods and datasets
To evaluate our method and compare it to other state-of-the-art works, we consider 34 functions gathered from the Nguyen (Uy et al, 2011) (noted N1 to N10), Keijzer (Keijzer, 2003) (K1-15), Vladislavleva (Vladislavleva et al, 2009) (V1-8), and Pagie (Pagie and Hogeweg, 1997) benchmark suites with varying levels of difficulty. The Nguyen benchmark suite is known to be the easiest one because it is mainly using one input feature and does not require optimizing for constant values in the expressions. Keijzer and Vladislavleva benchmarks are more complex as they contain functions with up to 5 inputs variables and also require to represent scaling constants. We also used Pagie (P1) (Pagie and Hogeweg, 1997) function, which has the reputation of being more challenging (McDermott et al, 2012). We applied several guidelines identified for Symbolic Regression benchmarking as provided by McDermott et al (2012). Our data generation procedure uses their functions and sampling intervals for train and test sets (McDermott et al, 2012).
The selected symbolic functions are benchmarked against the following grammar based methods: Grammar Guided Genetic Programming (G3P) (Whigham et al, 1995) is a grammar guided genetic algorithm build using the Deap library. (Fortin et al, 2012) Probabilistic Model Building Genetic Programming (GB-LGP) (Sotto and de Melo, 2017) updates the probability distribution of a grammar according to selected individual from an evolutionary population in gradient-descent like algorithm. As these methods are not directly available with our expression representation, we have re-implemented both methods. The following code is available on our Github repository.   (Sotto and de Melo, 2017). <e> is the start symbol, T = {<e>, <et>, <etw1>, ..., <varidx>} NT={x[], +,-,*,/,cos,sin, exp, log,1,..,nvar} To have comparable results between these methods, we propose to use the grammar described in Figure 4 for G3P, GB-LGP, and RBG2-SR (ours). It describes the transitions between 7 symbols and defines the action space to search into, made of at least 19+n var actions (with n var the number of features in the dataset). Datasets are generated using the drawing process detailed in Appendix A.
After hyperparameters search using a grid search method for all three algorithms, we choose the following best perfoming parameters for RBG2-SR: λ H = 0.005 and a learning rate α = 0.001. All methods are compared on a maximal horizon of 50 actions and each run is performed on a population/batch of 1000 expressions with a total of 2 millions expressions tested at most (corresponding to a 2000 iterations: batch size × nb training steps= 1000 × 2000 = 2M ).

Benchmark results
Expressions found by these methods on 30 independent runs are compared using Mean Squared Error (MSE) averages and standard deviations between the exact expression to uncover and best in-run solution of each algorithm. The last row of the table corresponds to the result of the Mann-Whitney U test for independent samples (Mann and Whitney, 1947). U-test results are summarised by counting the number of times each method performs better (symbol +), equivalently (∼), or worse (−) than other methods. "Equivalently" refers to the case where two (or more) methods are equally good, and have the same best results. A method is said to perform "worse" if at least one of the two other methods is performing "better" or "equivalently". Results for these benchmarks are shown in Table 1 with best results per function in bold and equivalent scores in italic. Table 1 Mean Squared Error and Standard Deviation scores for benchmarked methods, averaged over 30 runs (best results in bold). The symbol − is used when unable to compute a solution or when the solution error is larger than 10 10 . On the last two rows, is shown first the average MSE across all valid runs (out of 30 runs), and all benchmarks and then the count of times where each method is performing better (symbol +), equivalently (∼) or worse (−) than others performed using the Mann-Whitney U test    Fig. 5 Percentage of exact solutions found on the four benchmarks for the three methods RBG2-SR (blue), GB-LGP (orange), G3P (green). Better seen in color. 100% means that all 30 runs uncover the solution, 50% that 15 out of the 30 were able to match the right results As shown on the antepenultimate row in Table 1, our method lower the average MSE by one order of magnitude and performs statistically better than other methods on 22 out of 33 benchmarks according to the Mann-Whitney U test. The K5 benchmark is not taken into account in testing because all methods perform poorly and produce a high MSE error on this benchmark (above 10 14 ). We also note that we perform similarly to other methods for seven benchmarks, leading to an at least similar error in more than 90% of the tested benchmarks (30 out of 33 benchmarks).
More precisely, looking at the Nguyen benchmark, we show that our proposed method outperforms other methods for all functions, except for N4, N6, and N9, where the error is similar between RBG2-SR and G3P. Then, for Keijzer and Vladislavleva benchmarks, our method proposes solutions with a lower (resp. equivalent) error for 10 out of 15 (resp. 1/15) and 5 out of 8 (resp. 3/8) benchmarks. We also highlight that our method seems to complement the evolutionary G3P method, especially on the Keijzer benchmark, as G3P has a lower error on the few benchmarks where our method is weaker.
Let us now analyse the ability of our algorithm to precisely retrieve the exact target symbolic expression. It is also worth noting that our method reaches a zero error for several benchmarks (N1-3 and K7), meaning that we recover all the time the exact solution. More details on this result is described in Figure 5. In this Figure, we represent for each benchmark and method the number of times (in percentage) where we can recover the exact solution up to the numerical precision. First, we see that two of the three methods recover functions mostly on the Nguyen and Keijzer benchmarks. As we detailed above, unlike other methods for most benchmarks, the RBG2-SR method can often recover the exact solution for 9 out of 10 expressions of the Nguyen benchmark. On the same benchmark, G3P only recovers approximate solutions, and GB-LGP can only find four functions with a lower percentage. While GB-LGP recovers 3 out of 15 expressions from the Keijzer benchmark, our RBG2-SR method recovers 4 functions with a higher recovery rate. RGB2-SR also recover sometimes one function on the Vladislavleva benchmark. Note that these results are tied to the grammar used. The GB-LGP and G3P methods could potentially find more exact solutions by choosing another grammar. Moreover, among other unfound functions, some are difficult or even impossible to construct with the chosen grammar. This is especially true for functions such as K1-3 or V1 that require finding decimal constants, which is not supported in this version of the grammar.

Ablation study
In order to identify which elements of the RBG2-SR method are essential to the success of the expression search, we performed an ablation study on elements of the state definition s h and on the algorithm. In the algorithm itself, we try removing the risk-seeking objective J risk θ (keeping all trajectories) and the entropy loss term J entropy θ . Regarding the state definition, we compare state defined with and without: the current symbol σ h , the current mask m h , the current depth d h , the parent node a parent h , the siblings nodes a siblings h and the previously selected actions a past h . This ablation study uses the ten functions of the Nguyen benchmark and was run 10 times for each ablation/function combination. All results are compared to the baseline case where no element is occluded. The results are summarized in Table 2. First, when looking at the algorithm learning itself, we show in Table 2 that the combination of the two loss terms is relevant to our problem. The MSE significantly increases when removing one of these terms. The risk-seeking policy is crucial to this type of learning: removing the risk-seeking policy increases the error by 3200% when compared to the proposed method with risk-seeking. The entropy term is also of great importance, with an error increase of 860%.
We performed a second type of ablation on the state definition. In this part of the ablation study, we tried removing elements from state inputs from the Neural Network (the orange block called "state information"): either parent action information a parent , siblings action information a siblings , previously selected actions a past . Results from Table 2 indicate that these siblings, past actions, depth, and symbols are highly beneficial to the expression search since removing one of these terms increases the error by at least 63%. Depth and symbol information seems to be of almost equal importance for a successful expression search. Regarding parent information, it seems that removing this information tends to reduce the error by a small 4%. However, by looking more carefully at each benchmark, removing the parent information is only beneficial to the search for benchmark N9. Except for N1-3 and N7 (where both configurations always recover the exact expression), the complete proposed algorithm (Baseline) outperforms the parent ablation. For the remaining benchmarks, the corresponding increase lies between +50 and +55%. With these results, we choose to keep parent information in the state definition.

Experiment 2: Interpretability analysis of a use case
This work aims at describing an algorithm that provides interpretable symbolic solutions directly readable by humans. From the first experiment, we also see a potential application of our approach to more complex datasets with an unknown relationship between a set of observations X and a target variable y, where X and y variables may be of different physical units. In this scenario, the use of a grammar is particularly important to restrict outputs to realistic solutions in terms of physical units. For example it is not physically feasible to add a speed (measured in meters per second) with a distance (in meters).
Toward this physical interpretability goal, we designed a second experiment on the Airfoil Self-Noise dataset (Brooks et al, 1989), to compare the solutions proposed by our algorithm with the ones given by four state-of-the-art (Whigham et al, 1995;Sotto and de Melo, 2017;McConaghy, 2011;Petersen et al, 2021) methods. The dataset is accessible on the UCI Machine Learning Repository 2 . The tested methods are GB-LGP, G3P from the previous experiment, and two other non-grammar based methods: Fast Function Extraction (FFX) (McConaghy, 2011) which applies pathwise learning to a large set of nonlinear functions, and exploits the path structure to generate models that trade-off error/complexity. Deep Symbolic Regression (DSR from ) a deep learning algorithm which uses a RL approach to search the solution space without grammatical constraints.

Dataset Description
The problem we focus on in this experiment is predicting scaled sound pressure level (SSPL) on different size NACA 0012 airfoils. The estimation is made using the following variables: frequency (unit Hz), angle of attack (degree • ), chord length (meters m), free-stream velocity (meters per second m.s −1 ), suction side displacement thickness (meters m).The measurements are obtained using airfoils taken at various wind tunnel speeds and angles of attack. The span of the airfoil and the observer position is fixed for all measurements. The dataset is split between 70% in the train set and 30% in the test set.

Grammar construction
The first preprocessing step, before methods comparison, is the definition of a constrained grammar that contains premice of knowledge on the studied topic, such as the one used in Figure 6. The start symbol is <exp>. This symbol describes the dimensions (units) and structures the algorithm allows as a return. From the previous studies (Lau et al, 2009) on this data, we want to find an expression in the form: exp = constant − 10 * log10(child expression).
We also add several other structures to leave to the algorithm the freedom to explore.  To sum up the grammar from Figure 6, the three first lines describe what the units and non-units (composition of units) of the problem are. In the four final lines, the grammar constrains the operations on each dimension (or unit) to only physically consistent combinations. These lines also define how to go from one unit to another by using physical properties (such as velocity law). This part of the grammatical description is of particular importance to describe the expertise and knowledge we want to include to constrain the search space during the SR resolution.

Experiments and results
In this experiment, all algorithms are run 10 times (except for FFX (McConaghy, 2011) as it is deterministic), and the best expression of all runs is shown in Table 3. Their results are compared on the test set based on the MSE error, determination coefficient R 2 and complexity C. The complexity C is defined by the sum of operations and input features used in the formula. From Table 3, we can note that the best performing algorithm on this dataset is FFX, with an error of 10.3. However, this method also produces the expression with the highest complexity (around 140), above the maximal horizon H of 50. As seen in the column Expression, FFX solution is complex, not directly readable by a human and combine variables with different units. Among other four methods, we highlight that our proposed RBG2-SR method produces the second lowest error and highest determination coefficient of 0.72, while keeping

Conclusion and perspectives
This study proposes a new algorithm (RBG2-SR) for Grammar Guided Symbolic Regression using a Reinforcement Learning search approach that allows the inclusion of domain knowledge within the learning process and in the format of the solutions. We describe a POMDP modeling of the SR task in a grammatical action space. The proposed method is benchmarked against grammar-based state-of-the-art algorithms and shows significant improvements over other algorithms regarding the error metric and exact expression discovery. We performed an ablation study of the blocks of our algorithm and state definition. The results show that parent, sibling, past actions, depth and symbol informations are all important elements of the state definition. In the second experiment, we also show how the use of a grammar based approach could be useful and interpretable when working on a dataset with physical constraint between input features.
From the obtained results, we also foresee different perspectives to this work. First, when comparing G3P and RGB2-SR results, we envision improvements by doing cross-learning (Zhang and Zhou, 2021) between G3P to encourage exploration and our method for learning and sampling. Moreover, as the grammar construction process can be a time-consuming task, we could draw inspiration from techniques that automatically build ontologies (Emani et al, 2019) to automatically create and improve grammar. We also foresee application perspectives for our method to find interpretable policies that follow expert grammatical rules. From the second experiment, we also want to explore the behavior of our algorithm when dealing either with longer horizons up to 100 actions to create more expressive expressions or with shorter horizons for more concise and interpretable expressions. Code availability. The code is freely available on the Github repository whose link is shared in this document.