A Novel Neural Network-Based Symbolic Regression Method: Neuro-Encoded Expression Programming

TL;DR

This paper introduces Neuro-encoded Expression Programming (NEEP), a neural network-based symbolic regression method that creates a smooth, continuous search space for genetic programming, improving efficiency and accuracy.

Contribution

It proposes a novel neuro-encoding mechanism using RNNs to enhance genetic programming for symbolic regression, addressing issues of discrete solution space and landscape ruggedness.

Findings

Improved training efficiency on symbolic regression tasks.

Reduced test errors compared to traditional methods.

Smoother fitness landscapes facilitate better search performance.

Abstract

Neuro-encoded expression programming(NEEP) that aims to offer a novel continuous representation of combinatorial encoding for genetic programming methods is proposed in this paper. Genetic programming with linear representation uses nature-inspired operators (e.g., crossover, mutation) to tune expressions and finally search out the best explicit function to simulate data. The encoding mechanism is essential for genetic programmings to find a desirable solution efficiently. However, the linear representation methods manipulate the expression tree in discrete solution space, where a small change of the input can cause a large change of the output. The unsmooth landscapes destroy the local information and make difficulty in searching. The neuro-encoded expression programming constructs the gene string with recurrent neural network (RNN) and the weights of the network are optimized by powerful continuous evolutionary algorithms. The neural network mappings smoothen the sharp fitness landscape and provide rich neighborhood information to find the best expression. The experiments indicate that the novel approach improves training efficiency and reduces test errors on several well-known symbolic regression problems.

Tables

Table 1. Table 1: Test problems used in this paper. U ​ [ a , b , c ] 𝑈 𝑎 𝑏 𝑐 U[a,\ b,\ c] is c 𝑐 c samples uniformly randomized in [ a , b ] 𝑎 𝑏 [a,\ b] for the variable. E ​ [ a , b , c ] 𝐸 𝑎 𝑏 𝑐 E[a,\ b,\ c] are mesh points which are spaced equally with an interval of c 𝑐 c , from a 𝑎 a to b 𝑏 b inclusive.

Name	Variables	Function	Training Set	Testing Set
Sphere5	5	$x_{\_}{1}^2+x_{\_}{2}^2+x_{\_}{3}^2+x_{\_}{4}^2+x_{\_}{5}^2$	$U[1,11,1000]$	$U[1,11,1000]$
Dic1	10	$x_{\_}1+x_{\_}2+x_{\_}3+x_{\_}4+x_{\_}5$	$U[1,11,1000]$	$U[1,11,1000]$
Dic3	10	$x_{\_}1+\frac{x_{\_}2x_{\_}3}{x_{\_}4}+\frac{x_{\_}3x_{\_}4}{x_{\_}5}$	$U[1,11,1000]$	$U[1,11,1000]$
Dic4	10	$x_{\_}1x_{\_}2+x_{\_}2x_{\_}3+x_{\_}3x_{\_}4x_{\_}5+x_{\_}5x_{\_}6$	$U[1,11,1000]$	$U[1,11,1000]$
Dic5	10	$\sqrt{x_{\_}1}+\sin \left(x_{\_}2\right)+{\log }_{\_}e\left(x_{\_}3\right)$	$U[1,11,1000]$	$U[1,11,1000]$
Nico9	2	$x_{\_}{1}^4-x_{\_}{1}^3+x_{\_}{2}^2/2-x_{\_}2$	$U[-5,5,1000]$	$U[-5,5,1000]$
Nico14	6	$\left(x_{\_}5x_{\_}6\right)/\left(x_{\_}1/x_{\_}2x_{\_}3/x_{\_}4\right)$	$U[-5,5,1000]$	$U[-5,5,1000]$
Nico16	4	$32-3\frac{\tan \left(x_{\_}1\right)}{\tan \left(x_{\_}2\right)}\frac{\tan \left(x_{\_}3\right)}{\tan \left(x_{\_}4\right)}$	$U[-5,5,1000]$	$U[-5,5,1000]$
Nico20	10	$\sum _{\_}i=1^5\frac{1}{x_{\_}i}$	$U[-5,5,1000]$	$U[-5,5,1000]$
Poly10	10	$x_{\_}1x_{\_}2+x_{\_}2x_{\_}3+x_{\_}3x_{\_}4+x_{\_}4x_{\_}5+x_{\_}5x_{\_}6$ $+x_{\_}1x_{\_}7x_{\_}9+x_{\_}3x_{\_}6x_{\_}10$	$U[-1,1,250]$	$U[-1,1,250]$
Pagie1	2	$\frac{1}{1+x_{\_}{1}^{-4}}+\frac{1}{1+x_{\_}{2}^{-4}}$	$E[-5,5,0.4]$	$E[-4.95,5.05,0.4]$
Nguyen6	1	$\sin \left(x\right)+\sin \left(x+x^2\right)$	$U[-1,1,20]$	$U[-1,1,20]$
Nguyen7	1	$\ln \left(x+1\right)+\ln \left(x^2+1\right)$	$U[0,2,20]$	$U[0,2,20]$
Vlad3	2	$e^{-x}{x}^3\left(\cos x\sin x\right)\left(\cos x{\sin }^{2}x-1\right)\left(y-5\right)$	$x:E[0.05,10,0.1]$ $y:E[0.05,10.05,2]$	$x:E[-0.5,10.5,0.05]$ $y:E[-0.5,10.5,0.5]$
Energy	8	Energy efficiency of buildings
Concrete	8	Concrete compressive strength

Table 2. Table 2: Median, standard deviation and corresponding ranks of testing errors of the five compared algorithms. All differences are statistically significant according to a Wilcoxon test with a confidence level of 95 % percent 95 95\% . Symbols − - and + + represent that the proposed method is respectively significantly worse than and better than the other two methods (GP and GEP). The other cases are marked with = = .

	GEP	GP	GA-NEEP	PSO-NEEP	CMAES-NEEP
Sphere5	4.87e+04$\pm$6.58e+07	3.93e+02$\pm$1.40e+02	7.71e+02$\pm$4.46e+02	6.27e+02$\pm$4.35e+02	6.30e+02$\pm$1.17e+02
rank	5	1	4	2	3
			=	=	=
Dic1	6.00e+02$\pm$4.91e+07	1.55e+01$\pm$1.30e+01	2.04e+01$\pm$1.58e+01	2.83e+00$\pm$1.40e+01	4.97e-30$\pm$7.67e-02
rank	5	3	4	2	1
			=	+	+
Dic3	4.96e+02$\pm$5.58e+15	1.29e+02$\pm$2.61e+01	1.44e+02$\pm$2.39e+01	1.18e+02$\pm$2.51e+01	1.20e+02$\pm$1.30e+02
rank	5	3	4	1	2
			=	+	+
Dic4	7.00e+04$\pm$1.28e+11	6.68e+03$\pm$1.40e+04	3.01e+04$\pm$1.50e+04	5.12e+03$\pm$1.66e+04	3.62e+03$\pm$1.79e+03
rank	5	3	4	2	1
			=	=	+
Dic5	6.80e+00$\pm$1.35e+19	8.96e-01$\pm$8.70e-01	1.00e+00$\pm$2.61e-01	5.99e-01$\pm$2.50e-01	5.54e-01$\pm$1.21e-01
rank	5	3	4	2	1
			=	+	+
Nico9	4.80e+04$\pm$1.85e+19	3.64e+02$\pm$6.39e+03	1.32e+04$\pm$2.07e+05	1.27e+03$\pm$3.80e+03	2.65e+03$\pm$5.62e+04
rank	5	1	4	2	3
			=	=	=
Nico14	1.18e+07$\pm$1.47e+19	1.20e+07$\pm$6.95e+07	1.18e+07$\pm$1.12e+07	1.18e+07$\pm$1.93e+10	1.18e+07$\pm$5.89e+06
rank	3	5	1	4	2
			=	=	=
Nico16	4.47e+09$\pm$1.20e+18	4.48e+09$\pm$2.74e+11	4.46e+09$\pm$3.69e+11	4.91e+09$\pm$5.25e+12	4.46e+09$\pm$1.36e+13
rank	3	4	1	5	2
			=	-	=
Nico20	7.54e+02$\pm$1.21e+19	1.86e+03$\pm$4.64e+04	6.85e+02$\pm$2.15e+04	5.19e+02$\pm$5.52e+04	6.80e+02$\pm$1.54e+06
rank	4	5	3	1	2
			+	+	=
Poly10	5.48e-01$\pm$7.29e+00	3.24e-01$\pm$5.24e-02	3.21e-01$\pm$6.11e-02	3.21e-01$\pm$3.48e-02	3.17e-01$\pm$2.83e-02
rank	5	4	2	3	1
			=	+	+
Pagie1	9.61e-01$\pm$1.58e+19	1.26e-01$\pm$1.14e-01	1.95e-01$\pm$4.03e-02	1.24e-01$\pm$3.42e-02	1.21e-01$\pm$2.58e-02
rank	5	3	4	2	1
			=	=	=
Nguyen6	2.10e-01$\pm$2.57e+19	1.54e-01$\pm$1.61e-01	1.10e-01$\pm$1.22e-01	1.40e-02$\pm$3.28e-02	4.41e-03$\pm$1.47e-02
rank	5	4	3	2	1
			+	+	+
Nguyen7	2.63e-01$\pm$1.41e+19	3.92e-02$\pm$1.35e-01	3.30e-02$\pm$7.07e-01	2.19e-03$\pm$8.33e-02	1.15e-03$\pm$4.86e-03
rank	5	4	3	2	1
			+	+	+
Vlad3	7.63e+00$\pm$2.57e+20	1.22e+00$\pm$1.52e+12	9.47e-01$\pm$Inf	1.05e+00$\pm$Inf	1.01e+00$\pm$1.04e+33
rank	5	4	1	3	2
			+	=	+
Energy	1.06e+02$\pm$3.81e+18	2.58e+01$\pm$3.78e+01	4.52e+01$\pm$2.20e+01	2.17e+01$\pm$7.48e+00	2.34e+01$\pm$6.57e+00
rank	5	3	4	1	2
			=	+	+
Concrete	3.39e+02$\pm$5.78e+18	2.26e+02$\pm$6.63e+01	2.21e+02$\pm$2.66e+01	1.80e+02$\pm$3.41e+01	1.66e+02$\pm$3.38e+01
rank	5	4	3	2	1
			=	+	+
Avg. Rank	4.69	3.38	3.06	2.25	1.63

Equations

$$p_{\_}\mathrm{h}=\mathrm{round}(i_{\_}\mathrm{out}\cdot L+1),$$

$$p_{\_}\mathrm{t}=\mathrm{round}(i_{\_}\mathrm{out}\cdot (L-h+1)+h).$$

$$\{+,-,\ast ,/\}.$$

$$\{+,-,\ast ,/,\sin ,\cos ,e^n,\ln (\left|n\right|)\}.$$

Full text

††Aftab Anjum and Fengyang Sun —- Contribute equally to this article.

This work has been published on 2019 International Conference on Artificial Neural Networks, pp. 373–386, 2019. The final authenticated publication is available online at https://doi.org/10.1007/978-3-030-30484-331

1 Introduction

Symbolic regression (SR) [27] is to find an explicit function for simulation of user-defined data. Compared to implicit numerical (linear or nonlinear) regression analysis, SR can construct a function for a complex data without any prior knowledge and generally has powerful interpretability due to clear mathematic formula. Currently, genetic programming methods with linear representation [10, 3, 25] are mainly used to solve SR, which adopts a number of nature-inspired operators such as mutation and crossover to manipulate expressions and has presented decent performances on various applications. However, the linear representation approaches encode the expressions in discrete manner, which considers it as a combinatorial problem. Compared to continuous problems, the combinatorial problem do not provide sufficient and useful neighborhood information to aid searching [19, 28]. In addition, the local structure of the combinatorial problem are hard “sharp” style and the fitness landscape is not smooth where a minute change of the genotype can instill a substantial change of the phenotype [30, 8] and may cause oscillation of converging process [4] (Fig. 1). All the factors above make it difficult to find a desirable function fitting data for linear representation methods.

Recently neural networks have achieved great success in generative tasks [13, 23, 32]. In particular, neural networks have demonstrated considerable potential for generating texts or strings [2]. The genetic programming methods (such as gene expression programming) can decode a string to an expression tree that is equivalent to a mathematical function. Then the two facts make it possible to use neural networks to generate expressions, which converts the purely discrete encoding to continuous encoding and alleviates the aforementioned difficulties in solution space.

Therefore, we propose a neuro-encoded expression programming (NEEP), which constructs the mathematical functions with neural networks generating expression string. Instead of the discrete way, a small change in continuous weights vector only triggers similar form of function and makes slow-varying effects. Therefore, the continuous neural network mappings smoothen and soften the hard sharp discrete fitness landscape and provide more flexible local information. In this manner, the NEEP can adopt powerful continuous optimizers to finely adjust the weights of network and find better function.

2 Related Works

2.1 Symbolic Regression

The purpose of symbolic regression is to find an explicit function, which is the primary difference with numerical regression. For a predefined data, SR finds an explicit function $f:{\bf{x}}\to y$ that approximates the data with minimum error.

2.2 Genetic Programming

There are several types of methods that can be used to solve SR, e.g., analytic programming [39], fast function extraction [21], grammar evolution [26] and genetic programming [16]. Genetic programming methods (such as gene expression programming and standard genetic programming) are one type of the commonly used methods to solve symbolic regression. Standard genetic programming (GP) [16] tunes the tree structure of expression directly by nature-inspired operators, e.g., crossover is the exchange in subtrees of two chromosomes at certain nodes. GPs maintain good patterns but suffer from the explosion of tree size. Gene expression programming (GEP) [10, 36] constructs chromosomes with linear expression strings and provides an efficient way to encode syntactically correct expressions. For readers’ better understanding of our method, we introduce more details about GEP.

Gene expression programming encodes the expression tree structure into fixed length linear chromosomes. Structurally, GEP genes consist of head and tail. Head consists of function symbols and terminal symbols and tail contains terminal symbols only. GEP uses the population of linear expression strings, selects them according to their fitness values and introduces genetic variation through genetic operators.

The encoding design [40] has a significant influence on the performance of gene expression programming since it determines the search space as well as the mapping between genotypes and phenotypes. Traditional GEP adopts the K-expression representation [10], which converts a linear string to an expression tree by using a breadth-first travelling procedure. Li et al. [37]introduced new enhancements (P-GEP), which improve the encoding design by suggesting the depth-first technique of converting the string into the expression Tree. P-GEP increases the searching efficiency, but it is not scalable for complex problems. Automatically defined functions (ADF) [15] were, for the first time, introduced by Koza as a way of reusing code in genetic programming. Ferreira [9] introduced improvements (GEP-ADF) to encode the subfunctions into the expression tree which makes the GEP more flexible and robust. However, these encoding improvements are still based on discrete space, which lacks of sufficient neighborhood information. To the best of our knowledge, there are few approaches encoding the expression string in a continuous manner. Therefore, the difficulties from discrete encoding are yet to be solved.

2.3 Neural Network on Generation Task

Neural network has achieved success on generation task such as image [13], audio [23] and text [2] generation. In particular, several studies on text/string generation by neural networks are reviewed in this part. Bowman et al. [2] proposed an RNN-based variational autoencoder (VAE) language model that incorporates distributed latent representations of entire sentences in a continuous space and explicitly models holistic properties of sentences. Wang et al. [32] generated texts based on generative adversarial networks (GANs). This method builds discriminator with a convolutional neural network and constructs generator with RNN and VAE to solve the problem that GANs always emit the similar data.

To our knowledge, there is little research that uses neural network to solve symbolic regression by generating expression strings in spite of its success on generating texts. Liskowski et al. [18] proposed a method in which neural networks play a “pre-training” role to detect possible patterns in data, and then aids in finding the function rather than generating expression strings directly. Yin et al. [38] proposed a novel self-organizing reservoir computing methods by gene regulatory network. This method can process arbitrary sequences of inputs such as speech recognition, whereas generating expression is not a sequential problem and there is no external “formula” input for training. Another type of interesting works are word embedding methods [17], however in reality we do not have enough formula data to learn their underlying similarity.

3 Methodology

The neuro-encoded expression programming (NEEP) adopts an evolutionary algorithm (EA) to optimize the network connection weights, then use the neural network to generate the expression string, and the K-expression method to decode the string into a function and then calculate the simulation error (Fig. 2). The main contribution of this paper is the method for encoding and generating the expression tree, which is based on the output of a neural network.

3.1 General Framework

The first step of this method is initializing the parameters (e.g., pbest and gbest in particle swarm optimization [14], distribution mean and step size in covariance matrix adaptation evolution strategy [1]) and population (all the net weights to be optimized) of the evolutionary algorithm. When evaluating each network weights vector, the weights are inserted into a recurrent neural network that generates the expression strings composed of function operators (such as +, -, $*$, /) and terminal symbols (e.g., variables). After that, these strings are decoded into expression trees, which are equivalent to mathematical functions. By putting the data into the expression, we can compare its value and to the target value. We use the resulting error as the fitness value of each individual network in population. Then, we update all the necessary parameters in the evolutionary algorithm (e.g., update pbest, gbest and velocity according to fitness in PSO) and update the current solutions set (weights to be optimized). We repeat the above process until the termination condition is met (Fig. 2).

3.2 Encoding

Inspired by biological brain, neural network presents flexible learning and powerful representation capability [33] and has been widely designed and applied in various fields [34, 12]. The fully connected recurrent neural network [33] thus acts as the encoder which specifies the linear expression string. The neural network consists of hidden neurons and the output neurons, and all the hidden neurons are fully connected with each other. As is shown in Fig. 3, there are no formal input neurons in the network because no external information is input into the network during the generation of the string. Instead of back propagation, we use an evolutionary algorithm to optimize the weights between the output neurons and all the hidden neurons (see optimization part). The Gaussian shape function $f=e^{-x^{2}}$ is chosen as the activation function of the hidden and output neurons instead of a sigmoid function because of premature convergence during encoding the string.

Firstly, for symmetry, each neuron is initialized to zero. The hidden weights are uniformly randomized before the evolution and keep constant during the evolution. This fixed weights are the same for all the neural networks. As the behavior of the network can be chaotic, a small change in the initial condition could produce a significant difference in the later state [20]. As shown in [31], to reduce the instability we introduce fixed weights among all the hidden neurons and the remaining weights are still capable of finding the underlying pattern. All the other weights between the hidden and output neurons are uniformly randomized and are further optimized by evolutionary algorithms. After each time step, the output neurons are updated. Then, after all the time step the output neuron with the maximum value indicates which function or terminal symbol will be inserted into the expression string at a certain position. The position is determined by an additional neuron in the output layer in default, naming it the position insertion neuron. The process keeps inserting the symbols into the expression string until the desired length is achieved.

The string modeling is based on head and tail [10]. The number of output neurons is determined by the size of the function and terminal sets of specific problem. For position identification of head part symbols (terminals and functions), the position insertion in the head part of the string $p_{\_}\mathrm{h}$ can be assessed by Eq. 1.

[TABLE]

where $i_{\_}\mathrm{out}$ is output value of insertion neuron as position rate, which is distributed in (0, 1]. $L$ is the current length of the expression string. On the other hand, if $L$ is larger than the head size $h$ then, the corresponding terminal symbol will be inserted at a certain point of the tail part according to the value of position insertion neuron $i_{\_}\mathrm{out}$. The value of position insertion neuron in the tail part $p_{\_}\mathrm{t}$ can be calculated by the given Eq. 2.

[TABLE]

The whole process of symbol injection can be seen the encoding part in Fig. 2.

3.3 Decoding

The decoder is the translator which transfers the information from the string into the expression tree. The translation starting position is always the first position of the gene, whereas the last position of the gene does not necessarily coincide with the termination point.

Let us consider the encoded gene “\sqrt{\ }\+-**xxsinxyyyxyxxy” as represented in Fig. 2. This encoded gene can be translated into the expression tree by the breadth-first technique which is further decoded into the mathematical function. The fitness value of each mathematical function/expression is calculated by measuring how well it fits the data, using mean square error (MSE) between the predicted values and the desired values. The decoding process is the same as in GEP (see more details in [10]).

3.4 Optimization

In the NEEP framework, we can choose different evolutionary algorithms to optimize the neural network for producing the most accurate expression. Three versions of NEEP are proposed in this work, which is GA-NEEP (based on GA, genetic algorithm [11]), PSO-NEEP (based on PSO, particle swarm optimization [14]), CMAES-NEEP (based on CMA-ES, covariance matrix adaptation evolution strategy [1]). In all the evolutionary algorithms, the population (chromosomes or particles) are the weight vectors, their values are uniformly randomized and then insert into the neural network for encoding the expression strings.

We do not use back propagation (BP) because the problem is different from conventional supervised learning. The evolved function is not fixed during the calculation of derivative of weights. It is quite hard to obtain gradients with a uniform BP for all the problems and evolved functions. On the other hand, evolutionary algorithms [5] can conduct efficient optimization by approximating search direction without gradient information and become an appropriate choice of optimizer.

3.4.1 Genetic Algorithm

Genetic algorithm (GA) [11] is the search and optimization procedure inspired by the natural selection process of genetic. For evolving the population selection, crossover and mutation operators are used for creating better offspring than the parents. GA is a robust optimization method that is easy to implement and requires few problem information. In this study, the individuals are the weight vectors whose values are randomly created. Each weight vector is put into the corresponding neural network and then obtains its fitness value. GA is adopted in this work as a baseline method.

3.4.2 Partials swarm optimization

Partials swarm optimization (PSO) [14, 35] is a global numerical algorithm inspired by birds foraging. In PSO the particle migrate towards the direction of a combination of their former velocity personal best and global best and update its position. In this work, the particles are the weight vectors, and their values are randomly created. For the evaluation of these particles, they are inserted into the neural network and calculate their fitness values. PSO has shown robustness and efficiency in finding global optima.

3.4.3 Covariance matrix adaptation evolution strategy

Covariance matrix adaptation evolution strategy (CMA-ES) [1] is a leading stochastic and derivative-free method for continuous optimization of non-linear, non-convex functions. The CMA-ES samples candidate points with a multivariate normal distribution. It updates the mean vector and covariance matrix so that it encourages reproducing previously successful search steps based on the maximum-likelihood principle. The CMA-ES has shown advantageous convergence property when compared to many other evolutionary algorithms for a wide class of problems. CMSA-ES is adopted in this study as one of the optimizers because of the mentioned advantages and powerful performances.

4 Experiments

This section explores the performance of the proposed NEEP, and will not devise more sophisticated wrappers around GEP to improve the encoding way. The analysis will be limited to synthetic and benchmark regression problems. The three proposed methods are compared with standard GEP and standard GP. The problem configurations are outlined in Table 1, and the algorithm settings are described in the following subsection. Finally, the convergence behaviors and test accuracy are discussed.

4.1 Benchmark Configurations

We evaluated the proposed methods, GEP and GP on 14 synthetic benchmark problems [22, 7, 24] and 2 UCI data sets [6]. The Poly10 function is from [29]. All the benchmark problems are listed in Table 1. Function Poly10 and Sphere5 use the function set below

[TABLE]

The other functions use the function set below

[TABLE]

The division is protected by $f=x/(y+\varepsilon)$, where $\varepsilon$ is a very small number (e.g., 1E-100). Other benchmark details are listed in Table 1. All these difficult benchmark problems are commonly used due to their unique structural complexities with respect to objective formula. Several large scale benchmarks (e.g., 10 variables) for symbolic regression are considered one of the hard cases due to the difficulty of finding the solution in larger search space.

4.2 Compared Algorithm Configurations

Standard GEP and GP [16], which both are classic and powerful models in the field of genetic programming methods, are compared with the three instances of the proposed method (marked as GA-NEEP, PSO-NEEP, and CMAES-NEEP). For a fair comparison, all common parameters in the listed methods are initialized with the same value. All the algorithms in the experiments used a population size of 100, and the number of generations 500. Other parameters of GA, PSO and CMA-ES were specified by default. For GP, we used tournament size of 3, maximum tree depth of 10, maximum tree length of 61, maximum mutation depth of 4, maximum crossover of depth 10, maximum grow depth of 1 and minimum grow depth of 1. For GEP, we used header length of 30, a crossover rate of 0.7, mutation rate of 0.1, IS transposition of 0.1, RIS transposition of 0.1, and the inversion rate of 0.1. For the three proposed methods (GA-NEEP, PSO-NEEP and CMAES-NEEP), we used header length of 30, hidden neurons 40, time steps of 10, the initial fixed weights sparsity of 0.5 and the initial optimizing weight range of [-2, 2].

4.3 Results and Discussions

Table 2 summarizes the test errors obtained by GEP, GP and the three versions of NEEP on all the benchmark problems. The median and standard deviation are summarized over the 50 independent repeated trials for each of the 16 benchmarks function. It can be observed that the proposed methods (GA-NEEP, PSO-NEEP, CMAES-NEEP) significantly outperformed GEP and GP on 14 out of 16 problems according to the median of MSE, and perform competitively on the remaining problems. In particular, CMAES-NEEP reported dramatically lower MSE and more stable performance (according to their standard deviation values) on all the high dimensional data (Poly10, Dic1, Dic3, Dic4, Dic5, Nico20), while GEP and GP failed to locate the global optimum for these problems because the solution expressions of a high dimensional problem become overwhelming or extremely complicated. Therefore, such problems may become tough for the traditional GEP and GP due to their lack of capability to encode a complex function in a single string.

For the two regression data sets Concrete and Energy, the convergence curves in Fig. 4 and test errors in Table 2 reveal that CMAES-NEEP and PSO-NEEP have remarkable performance and high stability among all the compared methods. According to the convergence curves, in some functions (Nico16, Dic1, Dic3) these methods illustrate premature convergence and get stuck at a local optimum during evolution. For Nico9 and Sphere5 problem, GP sits as the best method among all the compared algorithms, due to its property of reusability of existing nodes during the encoding of the expression tree. On the other hand, GEP stands at the worst, and all other proposed methods show the competitive results concerning GP.

According to the test error results and the convergence curves, we can observe that the proposed method has shown more stable performance and faster training speed on most of the listed benchmark functions. In addition, we can obtain a rough conclusion that among the three evolutionary algorithms, CMA-ES is the most powerful, PSO is the second, and GA is the last one. These ranks also conform the general impression of their performances on many artificial benchmark functions in evolutionary computation. Therefore, it is important for NEEP to choose a strong optimizer for searching better neural networks.

5 Conclusion

This study proposes a novel continuous neural encoding approach to improve conventional linear representation in genetic programming methods for solving symbolic regression. Linear representation methods manipulate the expression tree structures in a discrete manner, which does not assist in a localized search of solution space. The neuro-encoded expression programming (NEEP) transforms the combinatorial problem to a continuous problem by using a neural network to generate an expression string, thus powerful numerical optimization method can be adopted to find a better mathematical function for symbolic regression. Empirical analysis demonstrates the method has the potential to deliver improved test accuracy and efficiency.

There are several interesting future research directions, such as to explore more neural network architectures for encoding and introduce the constant creation in string encoding mechanism. This new framework for now only improves one of linear representation methods and focuses on one application in spite of its potential for applying on more methods and applications. Therefore, one of the future works is to explore other types of genetic programming methods with neural networks. Another consideration is to apply NEEP to more applications (e.g., classification, digital circuit design and path planning).

Acknowledgments

This work was supported by National Natural Science Foundation of China under Grant No. 61573166, No. 61572230, No. 61872419, No. 61873324, No. 81671785, No. 61672262. Project of Shandong Province Higher Educational Science and Technology Program under Grant No. J16LN07. Shandong Provincial Natural Science Foundation No. ZR2019MF040, No. ZR2018LF005. Shandong Provincial Key R&D Program under Grant No. 2018GGX101048, No. 2016ZDJS01A12, No. 2016GGX101001, No. 2017CXZC1206. Taishan Scholar Project of Shandong Province, China.