### Genetic Programming Genetic Programming (GP) is the original implementation of Symbolic Regression (SR). Traditionally much of the literature has used the term Genetic Programming as both the problem and the implementation. Historically, term Symbolic Regression has been used as a sub-problem to Genetic Programming, as the task of evolving equations, and often as an evolutionary method for finding equations. We want to make a distinction between the task of symbolically regressing a parse tree in a grammar, from the method of realizing that tree. That is, separating the objective from the optimization method. From here on out, when we refer to Genetic Programming, we mean an evolutionary implementation of Symbolic Regression. And when we use Symbolic Regression, we mean the generalized regression problem, limiting our examples to equation regression, which we consider a specific sub-problem of Symbolic Regression. #### The Algorithm 1. **Overview** - Basics of Genetic Programming 1. **Representation** - Tree-based and beyond 1. **Evaluation** - Fitness measurement 1. **Optimization** - Genetic operators and selection 1. **Initialization & Termination** - It's important GP requires several components to be specified in its implementation. Representation and evaluation define the solution structure and how it is simulated to infer modeling accuracy. Initialization determines the search starting point and what portions of equation space are reachable. Selection and breeding policies determine which candidates are better and how the space of equations is searched. This section overviews each of these components in the GP algorithm. ##### **Overview of GP** The most common implementation of SR has been the evolutionarily inspired method called Genetic Programming (GP) [koza:1992:gen_prog](http://www.amazon.com/exec/obidos/ASIN/0262111705/geneticprogrammi). GP is part of the larger family of Genetic Algorithms (GA) a class of algorithms inspired by survival of the fittest' [ [holland:1962:outline](http://dl.acm.org/citation.cfm?id=321128), [goldberg:1988:genetic](http://link.springer.com/article/10.1023%2FA%3A1022602019183) ]. GP differs from GAs by solution representation. In GAs, solution representation is a fixed size structure. In GP, solutions vary in size, usually represented as a parse tree within a grammar. While GP can be applied to generic languages, and was originally intentioned to evolve computer programs, we shall restrict our discussion to mathematical equations. This is in line with the general focus of the GP literature. GP is a stochastic search heuristic over the syntax-tree representation of an equation. GP requires several components to be defined: representation, evaluation, selection, genetic operators, initialization & termination, and population management. Representation and evaluation, the first two components of a ML algorithm were detailed previously. We while expand upon their usage in GP here. The third component, optimization, is the combination of selection, the genetic operators, and population management. Initialization and termination are required because the success of GP is sensitive to their determination. The Basic GP Process begins with an initial, random population of models. Then each individual is evaluated to determine fitness. This is followed by a selection process for candidate survival and replication. Survivors are then recombined using methods which resemble sexual reproduction and mutation. GP continues this process for a number of generations, until a model of desired accuracy is discovered or a computational threshold is reached. For detailed overviews of tree-based Genetic Programming, see [ [koza:1992:gen_prog](http://www.amazon.com/exec/obidos/ASIN/0262111705/geneticprogrammi), [kouch:08:thesis](http://robotics.ee.uwa.edu.au/theses/2008-Genetic-Kouchakpour-PhD.pdf) ] For a nearly complete list of all GP publications, visit [GP-bib](http://www.cs.bham.ac.uk/~wbl/biblio/), maintained by Bill Langdon. ##### **Representation in Genetic Programming** Section [SR-components]("/sr/01b-sr/#components") described the various schemes for representing equations. GP uses the operators and operands of this framework to construct parse trees for equations. The parse trees are chromosomes of an individual. This analogy is important to the methodology which uses sexual reproduction, mutation, and natural selection to search. One note about representation, equation coefficients are real valued numbers. As we will see, having continuous parameters, and optimizing by evolution, makes the task much harder. **Representation Closure** In the basic GP algorithm, expressions are generated at random. This can result in invalid operations being evolved which violate our basic mathematical constraints. Closure is the property that an expression be valid mathematically. Invalid expressions include dividing by zero, taking the square root of a negative number, and violating other similar rules. A simple remedy exists in which invalid operations are removed and functions such as log are protected. This method has been termed *interval arithmetic* [ [keijzer:03:improving](), [kotanchek:2008:trustable]() ]. A second method is to place restrictions on the grammar, as in [ [hoai:2001:framework](), [hoai:2002:solving](), [hoai:2003:tree]() ]. Hoai prevents invalid expressions from being generated where the previous method removes them later. **Populations** are the pools of candidates equations GP has to work with. It is essential just a set of equations with nomenclature to reflect the biological evolution inspirations of GP. Basic GP uses a single population. Parents are selected from this population for breeding, children and parents compete for positions in the population between generations. **Schema Theory** developed as a framework for talking about partial solutions in GP [ [holland:1975:adaptation](http://psycnet.apa.org/index.cfm?fa=search.displayRecord&uid=1975-26618-000) ]. Originally associated with Genetic Algorithms, it has been useful in the discussions around solution propegation and convergence as well as the disruptiveness of crossover in GP. A schemata is a template for an individual where some of the parts are filled in and some are blanks'' and not yet determined. The Schema Theorem states that schemata which have fewer determined components and above-average fitness will increase exponentially in successive generations. This means that good schemata will propagate through the population. ##### **Evaluation in Genetic Programming** The basic GP algorithm uses the evaluation methods outlined in [SR-components]("/sr/01b-sr/#components"). The fitness of an individual is then determined by a combination of accuracy and size. We will revisit evaluation in GP later in the next section to discuss some variations on evaluation. ##### **Optimization in Genetic Programming** Optimization is performed over a series of iterations, know as generations in the GP scheme. Following the evolutionary motif, optimizing a population requires selection and replication. As the great MPUA says
The sole reason we exist is to survive and replicate. That's it. ~ Mystery
GP simulates the evolutionary process on equations. Replication uses *crossover* and *mutation* to create new individuals from existing ones. Selection uses the non-dominated sorting Pareto fronts described in Section [SR-???](). Here we will describe the particulars as the relate to GP. **Genetic operators** are the means by which GP searches the space of equations. The search uses crossover, an exploratory technique, to find unexplored regions and mutation, an exploitative technique, to search locally. A balance between both is needed to make progress and converge towards good solutions.
**Crossover** is analogous to sexual reproduction and uses two parent equations to produce child offspring equations. To perform crossover, subtrees of the parents are selected, snipped, and swapped. This method is similar to crossing DNA, except that GP is a branching representation, rather than linear. The original crossover method selects subtrees of the parents at random. Enhancements on selecting crossover points will be described below. It is considered a destructive operation due to the significant probability that good schemata of the tree are disrupted. Disrupted in this context means that a partial solution, stored in a subtree, has an internal node selected for the crossover point. To combat this, i.e. reduce the probability of disruption, candidate solutions grow extra 'genetic material' that has little value in the evaluation. This 'natural' growth, called bloat, is described below in the sections on limitations and enhancements. **Mutation** alters a single candidate equation to search locally around that candidate. Mutation simply selects a node of a tree and changes it to a different type. Examples of mutations are a change of variable from X to Y, an operation from addition to subtraction. Mutation is also the way the original GP algorithm optimizes real-valued coefficients. If a coefficient is selected for mutation, it scaled by a random value in some predefined range. **Injection** is a less used method in which a single candidate is altered by choosing a single point of change. It is similar to mutation, but instead of changing a node type, the subtree is trimmed and replaced by a freshly grown branch. In terms of exploratory and exploitative search qualities, injection lies between crossover and mutation. **Selection** is the most important part of evolutionary algorithms, and thus GP. Therefore it is very important to choose the selection mechanisms wisely [Dumitrescu:2000:EC](http://dl.acm.org/citation.cfm?id=350139). Selections role is to differentiate between candidates, based on their quality, for survival and replication. This artificial environmental pressure causes natural selection, raising the fitness of the population. [ [Dumitrescu:2000:EC](http://dl.acm.org/citation.cfm?id=350139), [goldberg:1991:comparative]() ] extensively cover the variations in selection schemes. Here, we will give an overview to show what is possible. There are two points where selection is used, for survival and replication. Replication is the chance of having genetics passed to an offspring. This is how good partial solutions are passed through the population and altered to find more of the solution. Parental selection can happen with or without replacement. Proportional selection favors better candidates by increasing the probability they will be selected. This can cause premature convergence to a local optima solution permeating the population, wiping out diversity. Rank based selection sought to overcome this by associating a position rather than a probability for chance of selection. Rank is often calculated by the dominate, dominated ration of a candidate. Rank selection benefits GP by increasing selection pressure when the population variance is low. Binary Tournament [[goldberg:1991:comparative]()] selects two individuals and competes them for one of the parental equations. This is down twice to generate both parents for a recombination. Binary tournament can be extended to compete N candidates against each other. Binary tournament is similar to selecting a pivot point in library implementations of quick-sort. Survival is the chance of an individual solution remaining in the population. In the elitist method, a parent is only removed if one of its offspring is better than it. In the pure method, all parents and children contend for a spot in the population. ##### **Initialization & Termination** **Initialization** determines the search starting point and what portions of equation space are reachable. GP uses random generation to create the initial population of candidate solutions. using a uniform random probability distribution There are two ways to grow a single parse tree. The first is to grow a full tree up to a specified depth. The second is to allow node selection to include variables. This causes generation to halt at unspecified times while traversing the tree, with the resultant tree having variable shape. To create the initial population, these two tree generation schemes can be used for the entire population or can be used for different percentages of the population. A method known as ramped half-and-half creates half the population with the full method and the second half with the variable grow method. A size proportional set can be created by growing an increasing number of trees as the maximum depth increases. For example, if we limit candidates to binary trees, we grow twice as many trees as the previous depth. Ensuring diversity of equations at the start is difficult. One method is two continually generate equations, eliminating duplicates and invalid expressions, until the population is filled. Another method is to run several short run GP searches, leveraging GP's randomness (non-determinism), to obtain a pool of partially optimized candidates. **Termination** Since Genetic Programming is a heuristic search, it has no intrinsic stopping point. Therefore, we must tell it when to stop, The most common method used is generational, i.e. upon reaching a predefined number of generations or iterations. This gives us a reasonable idea of how long a search will take as well as providing a measurement of what can be solved in a given number of steps. Another class is performance based stopping criteria. To use these, an accuracy metric needs to be defined. One such metric is simply the best performing candidates error. If the error is less than some defined threshold, then the search is terminated. This method, however, is sensitive to the scale of the data and error measurement. Another performance based method uses *hits*, the number of data points which a candidate equation is within some threshold of the real value [citation](). In this case, the search terminates if any candidate scores a hit on all data points or some percentage of all the data points. Other schemes attempt to measure the progress of the search. One such example is the number of generations since a new best solution has been found. If the search continues for a time without finding improvement, then termination occurs. No matter what means is used to determine the stopping point, the results of the search need to be returned. Generally a set of the best solutions are returned, providing some amount of diversity. This set usually comes from the Pareto frontier, but can also include solutions which come from successive frontiers. #### Limitations in GP 1. **Disruptiveness of Crossover** - the basic search operator has uniform selection probability. It is the reason for bloat and inefficiencies. 1. **Equation Bloat** - the growth in tree size to reduce probability of node selection. Bloat emerges to protect good partial solutions. 1. **Equation Redundancy** - results from multiple copies of the same equation. Compounded by commutative and associative properties. 1. **Loss of Population Diversity** - results when a partial solution dominates the population. Caused by bloat and good partial solutions. 1. **Premature Convergence** - results when a search converges upon a local optimum. Caused by loss of population diversity then next point. 1. **Fundamental Issus in GP** - topics which are inherient to the process or have been unaddressed by the GP community The basic GP algorithm has many limitations stemming from its evaluation, selection, and breeding sub-algorithms. Many research papers in the field have sought to address these limitations and offer solution. In this section we discuss these limitations, detailing where they come from. We relate the solutions which offer improvement, to each limitation, and discuss their effect on the GP algorithm. The limitations, as well as the counters, are ofter intertwined between several aspects of the GP algorithm. In fact, if the initial starting points do not give sufficient coverage of the search space, GP will not be able to reach areas of that space. GP has issues which it inherits from evolution, which is good at dealing with uncertainty and variability, but not the ideal method for honing in on the most fit individual. These issues stem from the difficulties in maintaining a diverse population of equations while converging towards the best equation. Representation, the genetic operators, selection mechanisms, and population structure are interwoven in complex ways. These difficulties, and advancements made in GP, have effects which permeate through the entire GP infrastructure and operation. Measuring and understanding these effects requires rigorous testing and analysis, something laking but being addressed [ [McDermott:2012:benchmarks](https://cs.gmu.edu/~sean/papers/gecco12benchmarks3.pdf) ]. ##### **Disruptiveness of Crossover** The basic version of crossover has uniform selection probability driven by random number generation. This is the main cause of bloat and inefficiencies. The diagram below shows the probability of selecting a node given the tree size.
Figure 4-1 - Probabilities on the Tree
### Fast Function eXtraction Fast Function eXtraction (FFX) is the first deterministic SR implementation [ [McConaghy:http:FFX](http://trent.st/ffx/), [McConaghy:2011:GPTP](http://trent.st/content/2011-GPTP-FFX-paper.pdf), [McConaghy:2011:CICC](http://trent.st/content/2011-CICC-FFX-paper.pdf) ]. FFX does not use genetic operators, tree based representations, or random number generation. Instead, FFX uses a Generalized Linear Model (GLM). and a set of basis functions derived from the input features. FFX then applies a series of regressions to fit the parameters of the basis functions to a desired complexity This makes FFX very fast, but also limits its flexibility and the complexity of solutions. #### Algorithm FFX uses a Generalized Linear Model (GLM) of the form: $$y = F(\vec{x},\vec{w}) = \sum \limits_b^B w_b f_b(\vec{x})$$ The GLM has linear coefficients to terms of the summation, and is a flexible version of ordinary linear regression. [ [McCullagh:1972:Paper](), [McCullagh:1989:Book]() ] The $$f_b(\vec{x})$$ are not required to be linear functions themselves, but rather linear in coefficients, to the terms of the summation. In other words, there are no coefficients inside any of the $$f_b(\vec{x})$$. To learn the coefficients of the GLM, FFX uses Pathwise Regularized Learning (PRL). PRL augments least squares by adding regularization terms and then sweeping across and returning multiple parameter vectors [ [ZouHastie:2005:Paper](), [Friedman:2010:Paper]() ]. PRL also has some interesting properties: 1. Learning speed is comparable or better than least squares 1. Can learn when there are fewer samples than coefficients 1. Can learn thousands or more coefficients 1. Sets of coefficient vectors are returned trading off the number of coefficients and training accuracy
Figure #: FFX Algorithm
{% highlight Python linenos %} def FFX(Data): bases = createBasisFunctions(features) for b in range(1,B): alpha = 1.0 while complexityOf(eqn) != b: eqn := PathwiseLearning(Data,alpha) adjustLambda(eqn,alpha) best.Push(eqn) return best {% endhighlight %} Following the psuedo code, FFX first creates a set of basis functions for the GLM. To do this, univariate bases from each variable with the operations ($$x^{\pm 0.5}, x^{\pm 1}, x^{\pm2}, abs(x),log(x)$$). In their example this produced 176 bases. Next, the univariate bases were combined to produce 3374 bivariate bases, resulting in 3550 total bases. By allowing bases to be in both the numerator and the denominator, the overall number of bases doubles to 7100. Then, for all $$b$$ from $$1 \rightarrow B$$, FFX derives a linear combination of $$b$$ basis functions. To learn a model, FFX applies pathwise regularized learning to fit the GLM coefficients. FFX will then adjust the $$alpha$$ value to change the emphasis in the Pathwise Learning. This has the effect of alternating between increasing and decreasing the number of bases, continuing until the number of function bases equals the desired model complexity. FFX repeats this process for a number of desired complexities so that a variation of models is returned. This enables similar behavior to GP returning the first Pareto Frontier. #### Limitations FFX works well for many problems, requiring far fewer computations GP, FFX, however, has its own limitations: (1) there are no coefficients or parameters within the bases, meaning more difficult, non-linear relationships are beyond its abilities. This issue could be addressed by using non-linear regression and abstract coefficients. (2) Equations such as $v = x*y*z$ are beyond the abilities of FFX. Individual terms of the summation are limited in complexity to a pair-wise combination of uni-variate and bi-variate bases determined at initialization. Seeding with increased basis functions could become prohibitive as the number of terms grows through pair-wise basis combinations. In the 13 variable example provided, the initial number of GLM basis functions was 7100. As FFX is incapable of finding many of the benchmarks, we did compare against it.
### Regression Techniques There are many methods for regression, which generally attempt to explain the relationship between variables. Here, we highlight several which are common or relavent to the subject matter. #### Linear Regression Linear regression is a technique for modelling a dependent variable from a set of independent variables and unknown parameters. The most common method for determining the parameters is least squares, which minimizes the squared residuals [ [Legrand:1805:orbites_des_comietes](https://books.google.com/books/about/Nouvelles_m%C3%A9thodes_pour_la_d%C3%A9terminati.html?id=FRcOAAAAQAAJ)]. The key is that the dependent variable is linear to the parameters of the basis terms constructed from the independent variables. Nonlinear functions map be used, so long as parameters remain linear. #### Nonlinear Regression Nonlinear regression is an extension to the linear variation where the independent variable can now have more complex, nonlinear relationships to the parameters and depenent variables. Unlike linear regression, no guarentee can be made to the optimality of the final parameter estimate. Most methods for estimating the parameters incorporate successive iterations which improve upon the accuracy. Signifcant algorithms include Gauss-Newton, Levenberg–Marquardt, and gradient descent [ [Fletcher:2000:Book](http://www.wiley.com/WileyCDA/WileyTitle/productCd-0471494631.html), [Madsen:2004:Methods](http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf)]. PGE uses the Levenberg-Marquardt (LM) algorithm is an iterative technique that finds a local minimum of a function that is expressed as the sum of squares of nonlinear functions. It has become a standard technique for nonlinear least-squares problems and can be thought of as a combination of steepest descent and the Gauss-Newton method. When the current solution is far from the correct one, the algorithm behaves like a steepest descent method: slow, but guaranteed to converge. When the current solution is close to the correct solution, it becomes a Gauss-Newton method \cite{lourakis04LM}. #### Support Vector Regression Support Vector Regression Machines (SVRM) ere proposed as a variation on the original Support Vector Machine (SVM) [ [Vapnik:1996:SVRM](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.10.4845) ]. In SVM, a classification boundary is sought which maximizes the margin between the classes. In SVRM, the basic idea is to flip around the SVM problem and find a line with where the deviation is minimized. SVM and SVRM work for both linear and nonlinear models and is a convex optimization problem [ [Smola:1996:Tutorial](http://lasa.epfl.ch/teaching/lectures/ML_Phd/Notes/nu-SVM-SVR.pdf) ]. #### Ridge, Lasso, Elastic Net **Ridge** regression is a regularization method for ordinary least squares. It uses the L2-norm of the parameter vector, imposing a penalization is added for excessively large parameters. Often a scaling term is used to control the effect of the penalty. [[Tikhonov:1943:stability](http://a-server.math.nsc.ru/IPP/BASE_WORK/tihon_en.html)] **Lasso** regression is similar to Ridge as a regularization method on ordinary least squares. Lasso uses the L1-norm, which causes terms to have zero valued parameters, thus effectively removing them from the model [ [Tibshirani:1996:lasso](http://www-stat.stanford.edu/~tibs/lasso/lasso.pdf), [stanford:http:lasso](http://statweb.stanford.edu/~tibs/lasso.html) ] **Elastic Net** is the linear combination of Ridge and Lasso techniques [[Zou:2005:elastic](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.124.4696)]. This method uses the offsetting features of both methods to find a sparse model with few non-zero terms.
### Graph Algorithms We include several graph problems and algorithms as they served as seeds of thought that grew into the PGE algorithm. #### Minimum Spanning Tree The Minimum Spanning Tree (MST) of a graph is a tree which minimizes the sum of edges which connect all vertices in a graph. MST is solvable in polynomial time using a greedy algorithm. **Prim's** algorithm is a greedy algorithm for finding the MST of an undirected graph. The algorithm starts with an arbitrary vertex and constructs the MST by adding the lowest weight edge of a node which is already in the MST. A priority queue is used to efficiently determine the next lowest cost edge to add. [ [Prim:1957:Shortest](http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6773228) ] **Kruskal's** algorithm is also a greedy algorithm for finding the MST of an undirected graph. The algorithm starts by treating each vertex as a separate MST. It then proceeds to select the minimum edge from all remaining edges. If the edge connects two distinct MSTs, then it is included, otherwise, the edge introduces a cycle and is subsequently discarded. A disjoint-set data structure can be used to effeciently determine this condition. [ [Kruskal:1956:Shortest](http://www.ams.org/journals/proc/1956-007-01/S0002-9939-1956-0078686-7/home.html) ] #### Single-source Shortest Path The single-source shortest path (SSSP) problem is to find the shortest path to all nodes, given a starting point. **Dijkstra's** algorithm for SSSP is a greedy algorithm which resembles Prim's. [ [Dijkstra:1959:algorithm](http://www-m3.ma.tum.de/twiki/pub/MN0506/WebHome/dijkstra.pdf) ]. At each step, the next closest vertex is added to the already visited set, thus expanding outwards in a greedy manner. A priority queue is used as well to improve upon runtime complexity. The difference is that at each step, the distances to vertices on the frontier are updated to reflect the latest addition, and thus possibly changing the order of the priority queue. Knuth generalized Dijkstra's algorithm to hyper-graphs [ [Knuth:1977:generalization](http://www.sciencedirect.com/science/article/pii/0020019077900023)]. **A\* search** is a graph traversal algorithm proposed as an extension to Dijksta's algorithm with better runtime performance. [ [Hart:1968:astar](http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=4082128)] The A\* algorithm cuts down on the size of the subgraph that must be explored, by using a heuristic to estimate a lower bound on the "distance" to the target. Similarly, it uses a priority queue to determine the next vertice to include. The difference is that the current distance is added to the estimated distance to the target. ">next   index

## Symbolic Regression

Symbolic Regression (SR) is the task of deriving mathematical formula from observational data without any fore-knowledge of the domain or problem. In essence, this is the scientific process performed by a computer. Hypotheses are formulated, tested against the observations, and compared for explanatory value. This process mirrors our human notions for mathematical discovery in science and engineering.

### The Problem

The overall goal of SR is to produce models which find a balance in the trade-off between accuracy and complexity. The best models will therefore be both simple and explanatory, elucidating the dynamics of the system under study. In this way, an expert is freed to think about the larger and more complex aspects of that system, gaining insights from the results of SR.

The domain of SR is to model a dataset with a mathematical formula or equation. SR makes no a priori assumptions about the form of a model. In SR, the search encompases the space of all equations. Compare this with more traditional methods for regression. Linear and non-linear regression assume a model and fit the model’s parameters to the data on hand. More sophisticated methods like the Fourier transform and neural networks can model almost any signal or function. However, they don’t elucidate an understanding of the interactions and behavior of that signal or function. SR is often considered a generalization of traditional methods for regression. A SR implementation will generate, tests, and validated a massive number of equations. Thus it is flexible, capable of modelling almost any signal, and also explanitory, producing models which are easily interpreted by humans.

Because SR models data with equations, it also finds a vast array of application domains. Simple relationships are modeled with explicit functions. Time variant relationships are modeled with differential equations. Several differential equations can be combined into a set, forming a dynamical system. Relationships with vary in both time and space are modeled with partial differential equations. Invariants, or conserved quantities can also be modeled with an equation and its partial derivatives. These are all incantations of ideas from calculus and are prominent throughout the sciences and engineering. From cell biology to ecology, physics to astronomy, chemistry to pharmaceuticals, engineering and finance; all fields rely on mathematical notions and theory to produce solid foundations upon which decisions can be made.

Looking forward, a SR technology will. be a tool of great assistance to domain experts. The ability to automatically discover analytical models from our increasing volumes and complexity of data is paramount to continued scientific progress. The equation and formula models are human interpretable, precisely because they are a human conception and language for talking about functional relationships found in the world. It is the relation between humans, math, and science; that SR holds a special place. As a tool for helping humans make sense in a data driven world. SR holds the potential to accelerate the pace of scientific breakthroughs.

#### The Search Space

The SR search space is the space of all equations. This search space is defined by a grammar through its terminals, non-terminals, and production rules. It is an exponential space that can be represented as a tree or a graph. As a tree, the space is the ordered expansion of applied productions, with the start symbol at the root. The graph is similar, but connects nodes in the tree that are equivalent. In the graph, there are multiple paths to an equation, coinciding with the multiple derivations for an equations parse tree.

Defining the building blocks from which equations can be built, also defines the space of representable expressions that be searched. The size of this space is infinite even when disregarding the massive increase brought on by real valued coefficients. This is because a grammar is infinite through its production rules.

Consider the following spaces:

1. The space of all polynomials in one variable
2. The space of all polynomials in more than one variable
3. The space of all polynomials with trigonometric functions
4. Increasing the depth of expressions (allowing parentheses ()’)
5. Allowing division (put any of the previous spaces on top, bottom, or both)

The complexity of the search space is also linked to how the search space is represented. Tree representations create a large fanout ratio in the search space, with repetition. Using dynamic programming, that repetition can be detected and the search space simplified. The next section will describe various representations used in SR. Subsequent chapters, especially the PGE Chapter, will discuss the search space and these issues further.

### The Components

Like all machine learning tasks, SR can be logically broken into three components: Representation, Evaluation, and Optimization [Pedro:20??:ml, Mustafa - Learning from Data]. Representations are the conceptual and implementation models we choose, and which define the hypothesis space. Evaluation determines a model’s fitness, and is defined by the objective function. Optimization is then the technique which explores the representational space and maximizes model fitness. Often these three components have a high degree of coupling or co-dependence. That is, optimization techniques usually have a representation and objective in mind.

In our work, we additionally include a memoization component for tracking the sub-problems previously solved. Memoization is necessary in SR as the representational space contains an extremely high degree of overlapping subproblems. If we did not detect overlap, we would experience a high degree of repetition and process the same equations many times. To our knowledge, this is a first for the SR research community.

#### Representation in Symbolic Regression

In SR, equations need to be represented in cilico. The choice of representation determines the methods we can later use to work with the equations. As we will see, there are tradeoffs when choosing between them. There is also opportunity to use multiple representations simultaneously and will prove to be beneficial.

Of the many equation representations which have been explored in SR research, binary trees are the most common. They fit naturally with the theoretical concepts and simplify the underlying algorithms. Tree representationss in SR can be subdivided into binary and n-ary. Binary trees are natural to computers, n-ary tries are natural to algebra. Graph and linear representations have been used as well. All of the representations are interchangeable and we are able to transform between them. This is important as each representation has its benefits and can be useful to different subtasks in the SR problem.

Linguistic Foundations

Before we explore each representation, some context for the discussion is in order. In this work, we take a strong point of view from a language perspective. This perspective mainly comes from a computational background, in terms of regular languages, but also from and algebraic perspective. It was this perspective which lead to the PGE algorithm and underpins much of our foundation.

From a programmers point of view, a grammar defines a language in a recursive manner [ Aho:1972:TPT, Dragon Book ]. A grammar is comprised of a start symbol, nonterminals, terminals, and production rules. The language generated by a grammar G, denoted L(G), is the set of sentences generated by G. A sentence S is in L(G) if there is a sequence productions which results S. This sequence of productions constructs the derivation for S. If a sentence S has a derivation then it is in the language G. Thus, grammar allows us to determine if a sentence is in a language.

Below is a simplified example for a grammar for mathematics.

START -> E
E -> E + T | E - T | T
T -> T * N | T / N | N
N -> Cos(E) | Sin(E) | Tan(E) | Log(E) | Exp(E) | L
L -> (E) | -(E) | (E)^(E) | TERM
TERM -> Constant | Variable


When working with languages in computers, expressions are represented as trees. The trees represent the application of the production rules for a grammar. They are often refered to as syntax or parse trees. The grammar can then be used to determine if a tree is valid by finding a matching set of production rule applications. In PGE, the production rules will also provide a means to recursively enumerate sentences in a language.

Binary Trees

In Symbolic Regression, equations are constructed as binary trees made of basic components called building blocks. The tree of building blocks makes up the DNA of an equation. It matches the parse or derivation tree for a given equation, and can also be called its Abstract Syntax Tree (AST). Building blocks come in two types, operators and operands. Operators are the internal nodes of the AST and include familiar mathematical functions $(+,-,*,/,sin,cos,...)$. Operands are the terminals or leafs of the AST and represent state-variables, constants, and real-valued numbers. There are two ways to represent constant, indexed and floating point. The real value floating point representation has been the most common in GP literature. The value is randomly assigned, mutated by multiplying it by a small number or swapping it out with another during crossover. Indexing associates a constant in the tree with an element into an array. This permits the coefficients to be optimized outside of the equation form. This can be done by a second evolutionary algorithm or by standard regression techniques. We use real-valued coefficients in GP and indexed coefficients in PGE.

Algebra Tries

The algebra trie, or just plain trie, is a tree structure where each operand can have a variable number of children. The algebra trie representation only affects two operators in SR, addition and multiplication. In fact, in any tree based representation, some operators (internal nodes) such as $cos$ and $log$ always have only one outgoing edge to their operand. Figure \ref{fig:eqn-trees} shows the comparable representations of several equations as both binary and algebra tries.

Figure # - Binary Trees and Algebra Tries

The algebra trie representation has no effect on the constituents or modeling ability of the equation. It does however, alter the equation tree’s size and structure. Changing size of equations effects the trade-off between parsimony and accuracy. This multi-objective trade-off, often called Pareto sorting, is the multi-objective optimization search of SR. Any change to the fitness function undoubtedly effects the operation of any SR implementation.

Additionally, the trie representation eases the design of simplification algorithms. These algorithms work within the tree to rewrite it to a canonical form. The simplification algorithms are analogous to rewrite rules from parser literature [ Aho:1972:TPT ], and are reminiscent of high school algebra.

The implications of tries and associated algorithms have largely been unexplored in the SR community. Despite the lack of detailed analyses, our theories and experience have shown their use to greately increase the tractability of SR. As such, we will use them exclusively and omit the use of binary trees.

Linear Forms

Linear representations come in a few flavors and are usually the result of a transformation. Linear GP considers that a program tree is compiled into a linear set of machine instructions. In Linear GP, the evolution occurs on the machine code by modifying the sequence of instructions [ Refs Here ]. In Fast Function eXtraction [ McConaghy:2011:FFX ], equations are represented as a linear combination of simple basis functions. FFX uses non-GP methods to determine which bases to use and what the coefficient should be.

Graph Representations

Graph representations encode the same information as tree representations. The difference is that a subtree which appears more than once in the tree, can have a single instance connected to by different parts of a graph representation [ Hod:2007:graphrep]. A second variation is to maintain a library of small trees. Then, indexing nodes reference which basis to use and where to use it [miller:2000:cartesian]. With graph representations it is possible to simultaneously update several parts of a model, to its benefit or detriment.

Restrictions and Simplifications

Simple restrictions disallow invalid mathematics, operations such as dividing by zero or taking the logarithm of a negative number. This is sometimes called interval arithmatic in the literature \cite{keijzer:03:improving}. In general, we wish to only consider valid equations.

Simplification groups like terms together, remove sub-expressions which evaluate to a constant value ($\sin\pi$, $\frac{x}{x}$, or $x-x$), and reduce unnecessarily complex expressions such as $\frac{1}{\frac{1}{x}}$. These algorithms, known as tree rewriting systems [ Joshi:1975:TAG, dershowitz:1982:orderings, huet:1980:equations ], convert one parse tree to an equivalent and simpler tree. There is debate as to how simplification effects the SR process [ kinzett:2008:using, kinzett:2009:online, mcree:2010:symbolic ]. Simplifications effect the structure of the AST through its manipulation, which in turn effects the areas of the overall search space. In addition, depth and size bounds are often placed on the equation’s derivation tree. Placing bounds on the tree explicitly places bounds on the search space. Generally, a small bound on the size or depth is started with and increased if the models returned are unsatisfactory.

#### Evaluation in Symbolic Regression

In SR, equations need to be evaluated to assess their accuracy and fitness. This involves processing on the training data, measuring the error, and then constructing a fitness value to be used for relative comparison.

Point evaluation is evaluating a model representation on a data point. This is generally a straight forward process and is repeated for every point in the training data to produce the predicted values.

Accuracy metrics are methods for estimating the residual error of a model. Error is estimated by comparing the predicted values with the true values using a standard measure. Typical measures include:

1. Mean Absolute Error (MAE)
2. Mean Squared Error (MSE)
3. Root Mean Absolute Error (RMAE)
4. Root Mean Squared Error (RMSE)

Less direct accuracy calculations include: total hits, which tallies the number of data points that an equation is ‘close’ to, determined by a threshold. Ranking calculates a score based on relative order by considering the Pareto dominates, is dominated ratio. The most common accuracy metrics are mean squared error and average absolute error.

Fitness values are used to determine the relative fitness of individuals as opposed to the absolute fitness. In SR, they are generally constructed from the values of competing objective functions. The most typical is to minimize complexity while minimizing error. Complexity is usually measured as the number of nodes that are contained in the tree.

The values used for fitness can further be manipulated prior to model comparison. The values can be weighted to express relative emphasis during optimization and skew the balance. Normalization of the values, on a per value basis, limits the effect of scale between fitness components. The next section describes how fitness value are used to make the determination of relative fitness.

#### Optimization in Symbolic Regression

Optimization is then the technique which explores the representational space and maximizes model fitness. This phase is highly related to the method for solving SR. Arguably, one could say this is the most significant point of difference between the various methods. The next section describes the three main methods for SR.

Figure # - SR Loop

Exploration operators

Exploration operators determine how the search moves around the search space. The method of exploration is tied to the method of representation, since the former must act upon the later. Generally, an exploration operator will input one or two models and output one or more models.

Exploration operators can be further subdivided into exploration and exploitation operators. Exploration operators make larger moves within the search space, while exploitation operators search locally around good partial solutions.

Selection Mechanisms

Selection in SR chooses a subset of models from a larger group or population. Many methods for selection exist, including tournament, fitness-proportional, and stochastic sampling [ Goldberg:1989:book ]. SR can also be a multi-objective optimization (MOO). In this case the Pareto objective can be used for selection [ fonseca:1993:genetic, horn:1994:niched ]. Selection is used to determine which models will undergo a subsequent process or which are privileged enough to pass on.

Pareto Selection

In multi-objective optimizations, such as SR, no single model metric is better than others. There is a trade-off between the various metrics. In SR, there are two metrics, a real-valued error and a discrete size. We prefer SR to find the simpler equations which model the data accurately.

The Pareto non-dominated sort, usually just called the Pareto front, addresses this trade-off between opposing objectives [ luke:2002:lexicographic, smits:2005:pareto, van:1998:evolutionary, fonseca:1993:genetic, horn:1994:niched ].

At the Pareto frontier, no one solution is better than the solutions in both complexity and performance. In mathematical terms, vector $\vec{u} \prec \vec{v}$ if $\forall i, f_i(u_i) \leq f_i(v_i)$. In the case of SR, we have first consider the discrete objective, parsimony, and then the continuous objective, accuracy. A smaller function will always dominate any function which is larger and less accurate. A large function $g$ dominates $f$ when its error is smaller. That is to say that we accept larger functions when they are more accurate and smaller functions are allowed to have more error.

Figure # - Pareto Fronts

#### Memoization

Memoization is a technique which records the results of expensive computations, and later returns the cached value when requested. A primary example is computing an iterated function such as the Fibonacci sequence. Memoization is useful in these recursive cases because the repetition of computation tends to grow as the recursion deepens. [ Wiki ] Dynamic Programming is a technique where a problem is broken down into smaller sub-problems which are solved independently and memoized. DP is available when there is an overlapping and optimal substructural patterns within the problem.

SR is such a problem and can make use of memoization. To date, we are unaware of any work beyond our own which has done so. In Chapter 4, we will see how the SR search space can be constructed to fit the DP paradigm. This will prove to be exceedingly beneficial.

### Implementations

Similarly, as in any machine learning task, there are several algorithms which try to solve the SR problem. Originally, SR was born from Genetic Programming (GP), an evolutionary algorithm for attempting this problem. The use of GP as a method for SR has been the de-facto until only recently. It is for this reason that SR as a problem is just starting to find itself researched outside of evolutionary algorithms. With the advent of Fast Function eXtraction (FFX) and Prioritized Grammar Enumeration (PGE), there are now competing methods to GP which are also deterministic.

SR implementations are the assemblages of algorithms and data structures for representing, generating, and testing hypotheses. The most common implementation of SR has been the evolutionarily inspired method called Genetic Programming (GP). More recently, Fast Function eXtraction (FFX) was proposed as the first deterministic SR implementation, based on ridge and lasso regressions. The main contribution of this work is Prioritized Grammar Enumeration (PGE), a second deterministic implementation, with foundations in dynamic programming. We will briefly describe each here with further details in subsequent chapters.

#### Genetic Programming

GP is a stochastic search heuristic over the syntax-tree representation of an equation. GP requires several components to be defined: representation, evaluation, selection, genetic operators, initialization & termination, and population management. he Basic GP Process begins with an initial, random population of models. Then each individual is evaluated to determine fitness. This is followed by a selection process for candidate survival and replication. Survivors are then recombined using methods which resemble sexual reproduction and mutation. GP continues this process for a number of generations, until a model of desired accuracy is discovered or a computational threshold is reached.

Figure # - GP Loop

#### Fast Function eXtraction

Fast Function eXtraction (FFX) is the first deterministic SR implementation [ McConaghy:http:FFX, McConaghy:2011:GPTP, McConaghy:2011:CICC ]. FFX does not use genetic operators, tree based representations, or random number generation. Instead, FFX uses a Generalized Linear Model (GLM). and a set of basis functions derived from the input features. FFX then applies Pathwise Regularized Learning to fit the parameters of the basis functions to a desired complexity This makes FFX very fast, but also limits its flexibility and the complexity of solutions.

Figure # - FFX Loop

#### Prioritized Grammar Enumeration

Prioritized Grammar Enumeration (PGE) is a deterministic algorithm for SR. PGE solves the SR problem from a language point of view, to prioritize the enumeration of expressions in that grammar. By recognizing the structure of the search space, PGE turns the SR problem into a graph traversal problem. PGE also recognizes multiplicity of points in in the search tree and overlaps them, resulting in the search graph. PGE uses dynamic programming techniques that memoizes previous results and prevent the compounding duplication of effort.

Initially, PGE starts with a set of basis functions. Each basis function is fit to the training data, evaluated on the testing data, and placed into the Pareto Priority Queue (PPQ). The PGE algorithm iteratively pops equations from the top of the queue for processing. Each of these equations is expanded by applying the grammar’s production rules recursively. The new equations are memoized using a specialized structure, the Integer Prefix Tree (IPT). The IPT uses an equation’s serialized representation to track which solutions have been processed already. The IPT ensures that we only consider an equation form once. This allows PGE to eliminate duplication of effort in evaluation and expansion. Unique equations are fit to the training data, evaluated on the test data, and pushed into the priority queue. PGE continues until a model of desired accuracy is discovered or a computational threshold is reached.

Figure # - PGE Loop

### Applications

Symbolic Regression can, and has been, applied to many problems from a multitude of domains. This is possible because the underlying tree or trie can represent multiple types of equations. The different types of equation can, in turn, model different problems. Here we describe the different types of equations and provide some sample applications.

#### Equation Types

Problem Type Target Equation Form
Explicit Equations $y = f(\vec{x})$
Differential Equations $\frac{dx_i}{dt} = f(\vec{x})$
Invariants $F(\vec{x}) = C$

Explicit Equations

Explicit functions are the most basic equation type. These are a direct mapping, or rule, from the independent input variables to the dependent output variable. Data is plugged in and the output is the result of depth-first evaluation. The most familiar of these is the line: $y = ax + b$. The line has one input variable and one output variable. Other familiar examples are the polynomial: $y = ax^2 + bx + c$ and the plane: $z = ax + by + c$. A real-world example is calculating housing price based on square footage and other features of a house.

Differential Equations and Dynamical Systems

Differential equations are functions which relate a quantity to its rate of change. They are prominent through the sciences. Dynamical systems are a set of differential equations, either heterogeneous or homogeneous. Each differential equations, in turn, is a rule for the rate of change, in one variable over time, to the current state of all variables. That is to say that the future depends upon the past.

To evaluate differential equations, techniques of numerical integration are required. There are many methods of numerical integration [kress:98:num_analysis]. Runge-Kutta 4 (RK4) is one such technique. It makes four smaller steps in time, updating the input variables along the way. Since RK4 involves four evaluations of a function, it requires four times the computational effort. RK4 also requires all of the differential equations, for a dynamical system, in order to simulate the system. RK4 integrates from current state to produce the next state of the system. It thus requires a temporal ordering of the data. We use the RK4 algorithm in this work to produce the data sets for the dynamical systems.

When performing an SR search, it is possible to decompose the differential equations for separate evaluation. We can approximate the integrated values of the other variables by substituting interpolated data. This partitions the SR search making the equation recovery simpler task [ hod:07:pnas ]. We call this Partitioned RK4 (PRK4). PRK4 requires the current and next states of the system, and four evaluations, in order to measure the fitness of an equation. It also maintains the temporal restriction of data in the evaluation process.

An alternative evaluation method for differential equations exists, if the numerical derivatives of the data can be calculated. The equations fitness is then measured by its ability to predict the numerical derivative. There are two reasons for doing this. The first is that small variations in the equation result in a larger magnitude of difference from the fitness function. The integrals of $x^2$ and $x^2+1$ diverge linearly, which means the integrated values of RK4 are very similar. The fitness landscape of the numerical derivative method is less smooth. The second reason for evaluating against the numerical derivative data is that we can perform a single point evaluation. We perform four times fewer evaluations and we can perform them at arbitrary points [ hod:08:mining ].

Invariants and Partial Differential Equations

The natural extension to differential equations is Partial Differential Equations (PDE). They are analogous in that they relate a function to rate of change. PDEs allow more than one independent variable, such as x, y, and z in three dimensional space. Invariant represent the solution to such a system system and are also called a conserved quantity. Well known examples from physics are the conservation of energy, mass, and angular momentum. Previous work has shown that it is actually possible to uncover the Hamiltonians for conservation of angular momentum, directly from measured data [ hod:09:science, hod:09:implicit_long ]. This begs the question:

Are there conserved quantities in biological or ecological systems.

Despite the potential for, PDEs and invariants are beyond the scope of this work.

top

### Further Considerations

In this work, we view SR as only a problem and view GP as an algorithm for performing SR. We attempt to make a distinction which the SR literature has not by using GP as a term for both the problem and solution. We view SR, not as a sub-class, but as a super-class to GP, as the regression within a grammar. A particular programming language is a specific grammar and an instance of the SR problem. So can be circuit or atomic regression. Mathematical equations are another instance of the SR problem and the one we focus on.

We want to reiterate the difference between problem and implementation. Symbolic Regression is the task which uncovers relationships in data by searching the expression space defined by grammar. An implementation is the assemblage of methods which realize the Symbolic Regression task. We believe the lack of separation has led to some missed opportunities is SR research. In this work, we step back from SR to fundamentally change the way we think about regressing symbolic expressions.

We also make this distinction because recent work [McConaghy:2011:FFX] introduces, to our knowledge, the first deterministic algorithm for SR called Fast Function eXtraction (FFX), offering a new means of performing an equation search.

#### Benchmark Problems

Benchmarking has been a problem for SR research. Candidate fitness metrics, methods for comparing implementations, and benchmark problems vary widely across the GP field. This situation has made effective comparison and progress difficult to assess. In order to compare SR implementations, and before trusting them on unsolved problems, methods for quantifying the effectiveness of their search capabilities are needed.

In 2012, McDermott:2012:benchmarks surveyed three years of literature from EuroGP and GECCO GP track, bringing this issue to the forefront of the community. Their aim was to start a discussion on unifying and standardizing the evaluation process in GP. We agree with these ideals and use 22 of their SR target functions for the evaluation of PGE. We do, however, disagree with the assumption in ~\cite{McDermott:2012:benchmarks}, that results should not be expected to be repeatable, and thus unverifiable by a third party. A non-GP practitioner will not likely use a tool which gives different answers each time it is used. This has been partially addressed by rate of convergence (how often an implementation finds an answer) and cumulative probability of success (the probability that an ideal solution would be found on or prior to generation $i$). Both of these methods require many trials. Nevertheless, we agree that on optimum is less obtainable and that a consensus needs to be reached on unbiased methods for comparison between different implementations.

#### Metrics for Success

Accuracy of Candidate Solutions

Accuracy predicts a candidates ability to model a given data set. Depending on the goals in using SR accuracy may or may not be the right measure. If we seek the most predictive model then accuracy can tell us which model is best. However, it is often the goal in SR to provide insight through meaningful relations. This is the trade-off between accuracy and parsimony. We want the simplest equation that isn’t to simple. In SR, we usually return multiple models instead of one single best model. This allows the domain expert can evaluate the varying degrees of complexity and accuracy.

It is usually the case that the most accurate model has also over fit the data set. Over fit models high frequency signals within the data by adding subexpressions to fit these signals. Validation sets overcome this issue by testing the model on unseen data after training. This results is an unbiased measurement, unless the learner is also being optimized. In that case, cross validation or a third data set is used as the never-before-seen testing data.

Rate of Convergence

The rate of convergence is defined as an implementation’s probability of finding the correct solution. This metric is because GP is a randomized algorithm which produces inconsistent results. A randomized algorithm’s results depend on the initial starting points. GP has non-deterministic behavior when creating the initial equation population and when genetic operations are applied during breeding. To compensate for the lack of consistency, GP researchers run many trials. The final results of all runs are tallied to come up with a probability for success. Rate of convergence can be a good measure when comparing the relative efficacy of GP implementations. It can also be used to tune the parameters of the GP implementation.

PGE doesn’t have a rate of convergence. It is a deterministic algorithm which executes the same way, in the same order, every time. Thus, given a set of parameters, PGE will either find the solution or it will not. The ability of PGE to find a solution depends upon the expansion functions and the amount of time alloted for searching.

Rate of convergence will not allow us to make direct comparison between GP and PGE. It will permit us to loosely infer problem complexity with the both implementations difficulty in finding the ‘correct’ solution. Further complicating this issue is that determining convergence is not well defined. The literature has used hit ratio, error thresholds, and $R^2$ correlation. Real convergence is the matching of form, which can be done by manual inspection or by tree comparisons. To our knowledge, no one has reported this type of convergence by directly stating how close the results are to the true form.

Is the solution present?

Determining if the exact solution is present only works for benchmark problems. Optimizing for this can mean the researcher has overfit their algorithm to the benchmark problems. Even so, an algorithm should return the answer to most benchmark problems.

The difficulties with this metric have been that GP is generally run multiple times to form the rate of convergence metric. It can become cumbersome to inspect the results for exact solution. Adding more difficulty to this process is the SR research has not considered the algebraic variations on an equation, that is, there is no canonical form. This makes the use of automated methods ineffective for detecting if the solution is present.

#### Metrics for Effort

Number of Iterations is often equated with how many steps it takes to find a solution. GP has a natural step between parents and children and is referred to as generations. Each generational step uses the current solutions and genetic operators to produce new candidate solutions. The number of generations is loosely tied to how long it takes a GP implementation to converge to an answer. In PGE, the generational step is applying the grammar productions to the best candidate solutions.

Measuring performance by iterations is easily a biased measure. It requires careful attention to detail. Different implementations use different logic and parameters, searching different amounts of equation space per iteration. Varying the number of equations or the layout of the population effects the search rate. Fewer iterations is not always better, but can be a useful measure to minimize when optimizing the parameters within a particular SR implementation. When comparing between implementations, iterations becomes misleading if other factors are not considered.

Clock Time

Clock time tells us how long, from when we hit enter, until we get our answer. All other considerations equal, it is an unbiased measure. Generally every thing is not equal, and why comparing several measures is useful and no measure should be used in isolation.

At some point, algorithms and computers become fast enough at a problem that time becomes less of an issue. To a domain expert, the difference between an hour and two, or hours and a day, may not be important. When they are modeling a large data set for a particular problem, a 24 hour turnaround is often acceptable, especially when one considers the time it takes to publish in a journal.

Time, despite being something not to fool with, is a powerful measure. Small incremental improvements compound in time to create meaningful progress. When time is coupled with other metrics, we can gain measures for the amount of work.

Number of Evaluations

Amount of work done is a good metric when we want the most return for what we put in. Work, however, can be difficult no quantify. In SR, work is usually measured by counting the number of evaluations made, and relating this measurement over time. Counting evaluations is easy in GP. Calculation is population size multiplied by generations and the number of data points. In PGE, the scheme is similar, however the number of models per iteration varies. A further unknown arises with the use of iterative fitting schemes like those used in FFX and with nonlinear regression. There is an extra iteration loop over the data points per model which needs to be accounted for.

Number of Models

Counting the number of candidates tells us how many equations we actually considered. If we count the number of evaluations we made on each equation, we get a measure of work. Most work occurs during the evaluation phase of the search process. The computation complexity of simulation and evaluation, far outweighs the time taken for the remainder of the search process. We can reduce this work by reducing either of the terms, number of equations or number of evaluations.

We can further refine this by considering the number of unique equations evaluated compared to the total number. We want to be smart about how many and which, equations and data points we choose. Coupling the unique to total equations with the number of evaluations will provide us a good comparison between the SR implementations. GP evaluates many equations fewer times, but experiences great amount of redundancy in form. PGE spends a larger amount of effort per equation, considers a form only once.

#### Reproducibility and the Importance of Determinism

This deterministic nature of FFX and PGE is an important feature for a SR algorithm. It is because of the link between human, math, and science. In science we expect reproducible results. We believe this is why GP has not become a predominant technology as predicted by its proponents, and likewise, SR for being so tightly coupled to GP for much of its history. Here, we take a strong stance that reproducibility is one of the key features to a SR implementation. Reproducibility provides the necessary reliability for scientific discovery and adoption of a new tool, as well as the foundation upon which research into SR can stand.

Deterministic execution has benefits. First, the algorithm is consistent and repeatable. This means that researchers are able to reproduce results easily. This is an important feature with large portions of research unable to be replicated [moraila:2013:measuring, MORE.. ]. Second, deterministic execution means that modifications to the algorithm can be studied in isolation with their effects becoming much easier to measure. Third, the PGE algorithm only needs to be run once on a given data set. This is a significant advantage over GP which requires 30 or more trials in order to obtain statistically significant results. In all of our modifications to the algorithm, we maintained PGE’s deterministic execution nature.

The following two quotes from prior publications are included in the hope that it will frame this work in the larger context of producing usable SR tools and technologies. Much of the research in SR aims at improving the GP implementation, with only limited results and inconsistent comparisons. By reconsidering the definition of the SR problem, we believe that PGE is an evolution in perspective and approach. Through determinism we may further our understanding of the problem and find ways to create such a true Symbolic Regression technology.

Genetic programming (GP) is not a field noted for the rigor of its benchmarking. Some of its benchmark problems are popular purely through historical contingency, and they can be criticized as too easy or as providing misleading information concerning real-world performance, but they persist largely because of inertia and the lack of good alternatives. Even where the problems themselves are impeccable, comparisons between studies are made more difficult by the lack of standardization. We argue that the definition of standard benchmarks is an essential step in the maturation of the field. We make several contributions towards this goal. We motivate the development of a benchmark suite and define its goals; we survey existing practice; we enumerate many candidate benchmarks; we report progress on reference implementations; and we set out a concrete plan for gathering feedback from the GP community that would, if adopted, lead to a standard set of benchmarks.
~ McDermott:2012:benchmarks

Outside the GP literature, SR is rare; there are only scattered references such as (Langley et al., 1987). In contrast, the GP literature has dozens of papers on SR every year; even the previous GPTP had seven papers involving SR (Riolo et al., 2010). In a sense, the home field of SR is GP. This means, of course, that when authors aim at SR, they start with GP, and look to modify GP to improve speed, scalability, reliability, interpretability, etc. The improvements are typically 2x to 10x, but fall short of performance that would makes SR a “technology” the way LS or linear programming is. We are aiming for SR as a technology. What if we did not constrain ourselves to using GP? To GP researchers, this may seem heretical at first glance. But if the aim is truly to improve SR, then this should pose no issue. And in fact, we argue that the GP literature is still an appropriate home for such work, because (a) GP authors doing SR deeply care about SR problems, and (b) as already mentioned, GP is where all the SR publications are. Of course, we can draw inspiration from GP literature, but also many other potentially-useful fields.
~ McConaghy:2011:FFX

### Genetic Programming Genetic Programming (GP) is the original implementation of Symbolic Regression (SR). Traditionally much of the literature has used the term Genetic Programming as both the problem and the implementation. Historically, term Symbolic Regression has been used as a sub-problem to Genetic Programming, as the task of evolving equations, and often as an evolutionary method for finding equations. We want to make a distinction between the task of symbolically regressing a parse tree in a grammar, from the method of realizing that tree. That is, separating the objective from the optimization method. From here on out, when we refer to Genetic Programming, we mean an evolutionary implementation of Symbolic Regression. And when we use Symbolic Regression, we mean the generalized regression problem, limiting our examples to equation regression, which we consider a specific sub-problem of Symbolic Regression. #### The Algorithm 1. **Overview** - Basics of Genetic Programming 1. **Representation** - Tree-based and beyond 1. **Evaluation** - Fitness measurement 1. **Optimization** - Genetic operators and selection 1. **Initialization & Termination** - It's important GP requires several components to be specified in its implementation. Representation and evaluation define the solution structure and how it is simulated to infer modeling accuracy. Initialization determines the search starting point and what portions of equation space are reachable. Selection and breeding policies determine which candidates are better and how the space of equations is searched. This section overviews each of these components in the GP algorithm. ##### **Overview of GP** The most common implementation of SR has been the evolutionarily inspired method called Genetic Programming (GP) [koza:1992:gen_prog](http://www.amazon.com/exec/obidos/ASIN/0262111705/geneticprogrammi). GP is part of the larger family of Genetic Algorithms (GA) a class of algorithms inspired by survival of the fittest' [ [holland:1962:outline](http://dl.acm.org/citation.cfm?id=321128), [goldberg:1988:genetic](http://link.springer.com/article/10.1023%2FA%3A1022602019183) ]. GP differs from GAs by solution representation. In GAs, solution representation is a fixed size structure. In GP, solutions vary in size, usually represented as a parse tree within a grammar. While GP can be applied to generic languages, and was originally intentioned to evolve computer programs, we shall restrict our discussion to mathematical equations. This is in line with the general focus of the GP literature. GP is a stochastic search heuristic over the syntax-tree representation of an equation. GP requires several components to be defined: representation, evaluation, selection, genetic operators, initialization & termination, and population management. Representation and evaluation, the first two components of a ML algorithm were detailed previously. We while expand upon their usage in GP here. The third component, optimization, is the combination of selection, the genetic operators, and population management. Initialization and termination are required because the success of GP is sensitive to their determination. The Basic GP Process begins with an initial, random population of models. Then each individual is evaluated to determine fitness. This is followed by a selection process for candidate survival and replication. Survivors are then recombined using methods which resemble sexual reproduction and mutation. GP continues this process for a number of generations, until a model of desired accuracy is discovered or a computational threshold is reached. For detailed overviews of tree-based Genetic Programming, see [ [koza:1992:gen_prog](http://www.amazon.com/exec/obidos/ASIN/0262111705/geneticprogrammi), [kouch:08:thesis](http://robotics.ee.uwa.edu.au/theses/2008-Genetic-Kouchakpour-PhD.pdf) ] For a nearly complete list of all GP publications, visit [GP-bib](http://www.cs.bham.ac.uk/~wbl/biblio/), maintained by Bill Langdon. ##### **Representation in Genetic Programming** Section [SR-components]("/sr/01b-sr/#components") described the various schemes for representing equations. GP uses the operators and operands of this framework to construct parse trees for equations. The parse trees are chromosomes of an individual. This analogy is important to the methodology which uses sexual reproduction, mutation, and natural selection to search. One note about representation, equation coefficients are real valued numbers. As we will see, having continuous parameters, and optimizing by evolution, makes the task much harder. **Representation Closure** In the basic GP algorithm, expressions are generated at random. This can result in invalid operations being evolved which violate our basic mathematical constraints. Closure is the property that an expression be valid mathematically. Invalid expressions include dividing by zero, taking the square root of a negative number, and violating other similar rules. A simple remedy exists in which invalid operations are removed and functions such as log are protected. This method has been termed *interval arithmetic* [ [keijzer:03:improving](), [kotanchek:2008:trustable]() ]. A second method is to place restrictions on the grammar, as in [ [hoai:2001:framework](), [hoai:2002:solving](), [hoai:2003:tree]() ]. Hoai prevents invalid expressions from being generated where the previous method removes them later. **Populations** are the pools of candidates equations GP has to work with. It is essential just a set of equations with nomenclature to reflect the biological evolution inspirations of GP. Basic GP uses a single population. Parents are selected from this population for breeding, children and parents compete for positions in the population between generations. **Schema Theory** developed as a framework for talking about partial solutions in GP [ [holland:1975:adaptation](http://psycnet.apa.org/index.cfm?fa=search.displayRecord&uid=1975-26618-000) ]. Originally associated with Genetic Algorithms, it has been useful in the discussions around solution propegation and convergence as well as the disruptiveness of crossover in GP. A schemata is a template for an individual where some of the parts are filled in and some are `blanks'' and not yet determined. The Schema Theorem states that schemata which have fewer determined components and above-average fitness will increase exponentially in successive generations. This means that good schemata will propagate through the population. ##### **Evaluation in Genetic Programming** The basic GP algorithm uses the evaluation methods outlined in [SR-components]("/sr/01b-sr/#components"). The fitness of an individual is then determined by a combination of accuracy and size. We will revisit evaluation in GP later in the next section to discuss some variations on evaluation. ##### **Optimization in Genetic Programming** Optimization is performed over a series of iterations, know as generations in the GP scheme. Following the evolutionary motif, optimizing a population requires selection and replication. As the great MPUA says
The sole reason we exist is to survive and replicate. That's it. ~ Mystery
GP simulates the evolutionary process on equations. Replication uses *crossover* and *mutation* to create new individuals from existing ones. Selection uses the non-dominated sorting Pareto fronts described in Section [SR-???](). Here we will describe the particulars as the relate to GP. **Genetic operators** are the means by which GP searches the space of equations. The search uses crossover, an exploratory technique, to find unexplored regions and mutation, an exploitative technique, to search locally. A balance between both is needed to make progress and converge towards good solutions.
**Crossover** is analogous to sexual reproduction and uses two parent equations to produce child offspring equations. To perform crossover, subtrees of the parents are selected, snipped, and swapped. This method is similar to crossing DNA, except that GP is a branching representation, rather than linear. The original crossover method selects subtrees of the parents at random. Enhancements on selecting crossover points will be described below. It is considered a destructive operation due to the significant probability that good schemata of the tree are disrupted. Disrupted in this context means that a partial solution, stored in a subtree, has an internal node selected for the crossover point. To combat this, i.e. reduce the probability of disruption, candidate solutions grow extra 'genetic material' that has little value in the evaluation. This 'natural' growth, called bloat, is described below in the sections on limitations and enhancements. **Mutation** alters a single candidate equation to search locally around that candidate. Mutation simply selects a node of a tree and changes it to a different type. Examples of mutations are a change of variable from X to Y, an operation from addition to subtraction. Mutation is also the way the original GP algorithm optimizes real-valued coefficients. If a coefficient is selected for mutation, it scaled by a random value in some predefined range. **Injection** is a less used method in which a single candidate is altered by choosing a single point of change. It is similar to mutation, but instead of changing a node type, the subtree is trimmed and replaced by a freshly grown branch. In terms of exploratory and exploitative search qualities, injection lies between crossover and mutation. **Selection** is the most important part of evolutionary algorithms, and thus GP. Therefore it is very important to choose the selection mechanisms wisely [Dumitrescu:2000:EC](http://dl.acm.org/citation.cfm?id=350139). Selections role is to differentiate between candidates, based on their quality, for survival and replication. This artificial environmental pressure causes natural selection, raising the fitness of the population. [ [Dumitrescu:2000:EC](http://dl.acm.org/citation.cfm?id=350139), [goldberg:1991:comparative]() ] extensively cover the variations in selection schemes. Here, we will give an overview to show what is possible. There are two points where selection is used, for survival and replication. Replication is the chance of having genetics passed to an offspring. This is how good partial solutions are passed through the population and altered to find more of the solution. Parental selection can happen with or without replacement. Proportional selection favors better candidates by increasing the probability they will be selected. This can cause premature convergence to a local optima solution permeating the population, wiping out diversity. Rank based selection sought to overcome this by associating a position rather than a probability for chance of selection. Rank is often calculated by the dominate, dominated ration of a candidate. Rank selection benefits GP by increasing selection pressure when the population variance is low. Binary Tournament [[goldberg:1991:comparative]()] selects two individuals and competes them for one of the parental equations. This is down twice to generate both parents for a recombination. Binary tournament can be extended to compete N candidates against each other. Binary tournament is similar to selecting a pivot point in library implementations of quick-sort. Survival is the chance of an individual solution remaining in the population. In the elitist method, a parent is only removed if one of its offspring is better than it. In the pure method, all parents and children contend for a spot in the population. ##### **Initialization & Termination** **Initialization** determines the search starting point and what portions of equation space are reachable. GP uses random generation to create the initial population of candidate solutions. using a uniform random probability distribution There are two ways to grow a single parse tree. The first is to grow a full tree up to a specified depth. The second is to allow node selection to include variables. This causes generation to halt at unspecified times while traversing the tree, with the resultant tree having variable shape. To create the initial population, these two tree generation schemes can be used for the entire population or can be used for different percentages of the population. A method known as ramped half-and-half creates half the population with the full method and the second half with the variable grow method. A size proportional set can be created by growing an increasing number of trees as the maximum depth increases. For example, if we limit candidates to binary trees, we grow twice as many trees as the previous depth. Ensuring diversity of equations at the start is difficult. One method is two continually generate equations, eliminating duplicates and invalid expressions, until the population is filled. Another method is to run several short run GP searches, leveraging GP's randomness (non-determinism), to obtain a pool of partially optimized candidates. **Termination** Since Genetic Programming is a heuristic search, it has no intrinsic stopping point. Therefore, we must tell it when to stop, The most common method used is generational, i.e. upon reaching a predefined number of generations or iterations. This gives us a reasonable idea of how long a search will take as well as providing a measurement of what can be solved in a given number of steps. Another class is performance based stopping criteria. To use these, an accuracy metric needs to be defined. One such metric is simply the best performing candidates error. If the error is less than some defined threshold, then the search is terminated. This method, however, is sensitive to the scale of the data and error measurement. Another performance based method uses *hits*, the number of data points which a candidate equation is within some threshold of the real value [citation](). In this case, the search terminates if any candidate scores a hit on all data points or some percentage of all the data points. Other schemes attempt to measure the progress of the search. One such example is the number of generations since a new best solution has been found. If the search continues for a time without finding improvement, then termination occurs. No matter what means is used to determine the stopping point, the results of the search need to be returned. Generally a set of the best solutions are returned, providing some amount of diversity. This set usually comes from the Pareto frontier, but can also include solutions which come from successive frontiers. #### Limitations in GP 1. **Disruptiveness of Crossover** - the basic search operator has uniform selection probability. It is the reason for bloat and inefficiencies. 1. **Equation Bloat** - the growth in tree size to reduce probability of node selection. Bloat emerges to protect good partial solutions. 1. **Equation Redundancy** - results from multiple copies of the same equation. Compounded by commutative and associative properties. 1. **Loss of Population Diversity** - results when a partial solution dominates the population. Caused by bloat and good partial solutions. 1. **Premature Convergence** - results when a search converges upon a local optimum. Caused by loss of population diversity then next point. 1. **Fundamental Issus in GP** - topics which are inherient to the process or have been unaddressed by the GP community The basic GP algorithm has many limitations stemming from its evaluation, selection, and breeding sub-algorithms. Many research papers in the field have sought to address these limitations and offer solution. In this section we discuss these limitations, detailing where they come from. We relate the solutions which offer improvement, to each limitation, and discuss their effect on the GP algorithm. The limitations, as well as the counters, are ofter intertwined between several aspects of the GP algorithm. In fact, if the initial starting points do not give sufficient coverage of the search space, GP will not be able to reach areas of that space. GP has issues which it inherits from evolution, which is good at dealing with uncertainty and variability, but not the ideal method for honing in on the most fit individual. These issues stem from the difficulties in maintaining a diverse population of equations while converging towards the best equation. Representation, the genetic operators, selection mechanisms, and population structure are interwoven in complex ways. These difficulties, and advancements made in GP, have effects which permeate through the entire GP infrastructure and operation. Measuring and understanding these effects requires rigorous testing and analysis, something laking but being addressed [ [McDermott:2012:benchmarks](https://cs.gmu.edu/~sean/papers/gecco12benchmarks3.pdf) ]. ##### **Disruptiveness of Crossover** The basic version of crossover has uniform selection probability driven by random number generation. This is the main cause of bloat and inefficiencies. The diagram below shows the probability of selecting a node given the tree size.
Figure 4-1 - Probabilities on the Tree
### Fast Function eXtraction Fast Function eXtraction (FFX) is the first deterministic SR implementation [ [McConaghy:http:FFX](http://trent.st/ffx/), [McConaghy:2011:GPTP](http://trent.st/content/2011-GPTP-FFX-paper.pdf), [McConaghy:2011:CICC](http://trent.st/content/2011-CICC-FFX-paper.pdf) ]. FFX does not use genetic operators, tree based representations, or random number generation. Instead, FFX uses a Generalized Linear Model (GLM). and a set of basis functions derived from the input features. FFX then applies a series of regressions to fit the parameters of the basis functions to a desired complexity This makes FFX very fast, but also limits its flexibility and the complexity of solutions. #### Algorithm FFX uses a Generalized Linear Model (GLM) of the form: $$y = F(\vec{x},\vec{w}) = \sum \limits_b^B w_b f_b(\vec{x})$$ The GLM has linear coefficients to terms of the summation, and is a flexible version of ordinary linear regression. [ [McCullagh:1972:Paper](), [McCullagh:1989:Book]() ] The $$f_b(\vec{x})$$ are not required to be linear functions themselves, but rather linear in coefficients, to the terms of the summation. In other words, there are no coefficients inside any of the $$f_b(\vec{x})$$. To learn the coefficients of the GLM, FFX uses Pathwise Regularized Learning (PRL). PRL augments least squares by adding regularization terms and then sweeping across and returning multiple parameter vectors [ [ZouHastie:2005:Paper](), [Friedman:2010:Paper]() ]. PRL also has some interesting properties: 1. Learning speed is comparable or better than least squares 1. Can learn when there are fewer samples than coefficients 1. Can learn thousands or more coefficients 1. Sets of coefficient vectors are returned trading off the number of coefficients and training accuracy
{% highlight Python linenos %} def FFX(Data): bases = createBasisFunctions(features) for b in range(1,B): alpha = 1.0 while complexityOf(eqn) != b: eqn := PathwiseLearning(Data,alpha) adjustLambda(eqn,alpha) best.Push(eqn) return best {% endhighlight %} Following the psuedo code, FFX first creates a set of basis functions for the GLM. To do this, univariate bases from each variable with the operations ($$x^{\pm 0.5}, x^{\pm 1}, x^{\pm2}, abs(x),log(x)$$). In their example this produced 176 bases. Next, the univariate bases were combined to produce 3374 bivariate bases, resulting in 3550 total bases. By allowing bases to be in both the numerator and the denominator, the overall number of bases doubles to 7100. Then, for all $$b$$ from $$1 \rightarrow B$$, FFX derives a linear combination of $$b$$ basis functions. To learn a model, FFX applies pathwise regularized learning to fit the GLM coefficients. FFX will then adjust the $$alpha$$ value to change the emphasis in the Pathwise Learning. This has the effect of alternating between increasing and decreasing the number of bases, continuing until the number of function bases equals the desired model complexity. FFX repeats this process for a number of desired complexities so that a variation of models is returned. This enables similar behavior to GP returning the first Pareto Frontier. #### Limitations FFX works well for many problems, requiring far fewer computations GP, FFX, however, has its own limitations: (1) there are no coefficients or parameters within the bases, meaning more difficult, non-linear relationships are beyond its abilities. This issue could be addressed by using non-linear regression and abstract coefficients. (2) Equations such as $v = x*y*z$ are beyond the abilities of FFX. Individual terms of the summation are limited in complexity to a pair-wise combination of uni-variate and bi-variate bases determined at initialization. Seeding with increased basis functions could become prohibitive as the number of terms grows through pair-wise basis combinations. In the 13 variable example provided, the initial number of GLM basis functions was 7100. As FFX is incapable of finding many of the benchmarks, we did compare against it.