Evaluation of Population Initialization Methods for Genetic Programming-based Symbolic Regression Lukas Kammerer1[0000−0001−8236−4294] , Gabriel Kronberger , Deaglan J. Bartlett2[0000−0001−9426−7723] , 3[0000−0003−0685−9791] Harry Desmond , Pedro G. Ferreira2[0000−0002−3021−2851] , and Stephan Winkler1[0000−0002−5196−4294]
arXiv:2606.31990v1 [cs.NE] 30 Jun 2026
1[0000−0002−3012−3189]
Heuristic and Evolutionary Algorithms Laboratory, University of Applied Sciences Upper Austria, Hagenberg, Austria 2 Astrophysics, University of Oxford, Oxford, UK 3 Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth, UK [email protected] 1
Abstract. We analyze the effect of optimizing the initial population of genetic programming (GP) for symbolic regression (SR) on the accuracy and complexity of solutions. We compare three well-established random initialization methods as well as initialization with small optimized solutions from exhaustive symbolic regression (ESR) using a GP/SR implementation which is based on the multi-objective evolutionary algorithm NSGA-II. We compare the final Pareto fronts found with each initialization method on twelve synthetic problems of varying complexity and one real-world dataset. We find no significant differences in accuracy or model complexity among the initialization methods. The initial advantage of initialization with ESR disappears after only a few generations. Our results show that, given similar diversity in the initial population, the effect of the initialization method in GP-based symbolic regression on the final Pareto front is negligible. Keywords: Symbolic regression · Genetic Programming · Population initialization · Tree creators.
1
Introduction
Symbolic regression (SR) is a machine learning problem where the goal is to find a human-interpretable mathematical expression of any possible functional form that best fits a given dataset [11, 12]. The most common approach to solve SR problems is genetic programming (GP), which evolves a population of individuals in the form of expression trees representing mathematical expressions. GP typically initializes a population randomly and then evolves it using selection, crossover, and mutation operations. This random initialization ensures a diverse set of individuals that can be used by the crossover operator to create new and potentially better models [11].
2
L. Kammerer et al.
Most research on improving GP/SR towards more accurate or shorter models has focused on the evolutionary process itself, for example by introducing new selection, crossover, and mutation operators [14]. Methods for creating an initial population and therefore initializing the starting condition of the evolutionary process are referred to as population seeding. These aim to provide aspects in the initial population that are beneficial for the evolutionary process, such as diversity [21, 10]. Given the crucial role of population diversity for the success of GP, most developments on population initialization focus on creating a uniform distribution of symbols or number of nodes in expression trees, such as the socalled grow method, the full method, the combination of the latter called the ramped half-and-half method [11], the probabilistic tree creators PTC-1 and PTC-2 [17] or the balanced tree creator (BTC) [6]. Different approaches to population seeding exist to deliberately bias initial populations. Methods that optimize the initial population’s fitness showed clear improvements in combinatorial optimization [21, 8, 10], while related attempts in GP for program synthesis [1] did not provide significant improvements over random initialization. Another approach by Langdon and Nordin [15] took overfit individuals as the initial population and used GP with parsimony pressure to identify well-generalizing models. Mundhenk et al. [19] used a recurrent neural network (RNN) to generate an initial population. After short GP runs that were seeded with the RNN, the best final models were then used to further train the RNN and generate a new population. This work analyses the effects of optimizing the initial population versus wellestablished random initialization methods. Our main hypothesis is that an initial population that is already strong in accuracy and complexity either accelerates convergence or improves the final Pareto front, provided that it retains sufficient diversity. To check this, we use short, well-fitting expressions produced by exhaustive symbolic regression (ESR) [3] for the initial GP population instead of random seeding. ESR performs a brute-force search over the space of all possible algebraically unique functional forms up to a certain number of nodes in the expression tree, which is referred to as complexity. Parameters are optimized numerically, and we assume in our experiments that the best parameters are found for each functional form. Due to its exhaustive nature, ESR is limited to a much smaller search space than GP/SR and covers only short, univariate models of complexity up to ten. We use the GP/SR implementation Operon4 [6] which has shown strong performance on SR benchmarks [14]. As baselines for comparison, we use the grow initialization method, PTC-2, and BTC. The grow method [11] creates trees of any complexity and shape within a given depth limit, which is the only available parameter to control the shape of the resulting trees. Due to this lack of fine-grained model complexity control, the grow method might generate only small or imbalanced trees that lead to low diversity, or too large trees [17, 4]. PTC-2 [17] creates trees from a customizable distribution of model complexity and operator occurrence. By default, Operon 4
specifically the Python wrapper PyOperon, version 0.6.0.
Evaluation of Population Initialization Methods
3
uses the balanced tree creator, BTC, which tends to build balanced trees of minimal depth and uniformly distributed complexities [6]. The use of ESR for population seeding is motivated by its guarantee to find the most accurate and algebraically unique models in their shortest form within the given complexity limits [3]. The population returned by ESR can be considered as an optimal population with respect to accuracy and complexity within Operon’s initialization complexity bounds because Operon also uses a complexity limit of ten for random initialization. ESR identifies and removes algebraically duplicate models in its search and returns a set of mathematically unique models, which leads to syntactically diverse models, as we demonstrate below. We compare ESR initialization with three different random initialization methods on twelve synthetic problems of varying complexity, and one real-world problem. Our results show that the initial advantage of the ESR-initialized population quickly diminishes during the evolutionary process within only a few generations. Therefore, there is no significant improvement in terms of accuracy and complexity. Only for simple problems, in which the ground truth is very close to the best model found by ESR, do we find a significant improvement with an ESR-initialized population over the random initial population. The experiments also show that all three random initialization methods perform equally well, which implies that the specific choice of random initialization method barely has an effect on the final Pareto front of models in GP/SR. Section 2 describes the experimental setup, including the initialization methods, the GP settings and the benchmark problems. Section 3 first outlines preliminary results regarding diversity and accuracy in the initial population to ensure a proper experimental setup. It then shows the modeling results of the synthetic problems and verifies their implications on a real-world dataset. Section 4 discusses and concludes our work.
2
Methodology
ESR is run before Operon to use the resulting set of best models as the initial population for Operon. For each dataset we run the ESR fitting phase with the same precomputed set of expressions from [2] that uses arithmetic operators and the power function. To match Operon’s distribution of model complexity in the initial population, we run ESR with different maximum complexities. As we use a population size of 1000 in GP and only binary operators that result in odd complexities, we run ESR with maximum complexities of three, five, seven, and nine and take the 250 most accurate models from each run. We use duplicate models for complexity three and five as for these complexity values fewer than 250 unique models exist. Given the deterministic nature of ESR (assuming perfect parameter optimization), we run it only once per problem. Since GP is a stochastic method, we perform 1000 repetitions of Operon with different random seeds for both ESR initialization and each random initialization method, as outlined in Fig-
4
L. Kammerer et al.
ure 1. We note that the ESR optimization process is performed in addition to GP and therefore uses additional computational budget compared to GP with computationally cheap random initialization. We compare the models from all resulting 1000 Pareto fronts of each method and plot the distribution of the error of the Pareto-optimal models from all runs over their respective complexity for all initialization methods. By comparing all distributions, we show whether changing the initial population leads to improvements regarding both accuracy and complexity in the final Pareto fronts. Since we are interested in the Pareto front of models over accuracy and complexity, we use the multi-objective NSGA-II [7] implementation of Operon. As objectives, we use the normalized mean squared error (NMSE) and the complexity, which is also referred to as the model length. The NMSE is defined as the mean squared error divided by the variance of the target variable in the training data. We use the same hyperparameters for all experiments, which are shown in Table 1, while only adapting the maximum complexity to specific problems. Problems 1 to 5 use a maximum complexity of 20, while Problems 6 to 12 use a maximum complexity of 50. The complexity limit for models in the initial population is set to ten, except for the grow method, which only supports a depth limit. Parameter values in randomly initialized models are sampled from a Gaussian distribution with zero mean and σ = 1. In ESR-initialized populations, we retain the parameter values that were found by ESR in each model. During the evolutionary process, parameter values in each model are optimized with the Levenberg-Marquardt algorithm [16, 18] for a maximum of 15 iterations. The function set is based on the functions used by ESR, with the only difference being that univariate inverse operations are converted into binary division operations with one as the numerator. The function set consists of arithmetic operators and the composite powerabs operator, defined as powerabs(a, b) = |a|b , which extends the power function to negative bases and ensures real-valued outputs for non-integer exponents. Operon’s default linear and variable scaling was disabled to identify functions in their shortest representation and to reuse ESR’s
for BTC/PTC-2/grow 1000 repetitions
Problem
Random Initialization
ESR-Initialization
GP
GP
1000 Pareto Fronts
1000 Pareto Fronts
1000 repetitions
Compare Distributions Fig. 1. We compare all initialization methods by running 1000 repetitions of GP with different random seeds on twelve synthetic and one real-world problem. We then compare the resulting Pareto fronts for each problem.
Evaluation of Population Initialization Methods
5
Table 1. GP/SR Hyperparameters used for all experiments. The maximum length depends on the problem. Other hyperparameters are set to their default values. Hyperparameter Value Algorithm: NSGA-II Objectives: NMSE, Length Population Size: 1000 Generations: 500 Selection: Tournament (size=3) Mutation Rate: 25% Max. Length: 20 or 50 Max. Initialization Length: 10 (only ESR, PTC-2 and BTC) Max. Initialization Depth: 5 Linear Scaling: off Variable Scaling: off Parameter Opt.: LM, max. 15 iterations Function Set: +, −, ×, ÷, powerabs Terminal Set: x, parameters
models in their original form. The most relevant parameters are shown in Table 1. All other parameters are set to their default values. We use twelve synthetic problems with varying complexity and a small noise level for comparison, as shown in Table 2. Each problem has a training set of size 100 and a test set of size 1000. We choose univariate synthetic benchmarks with known noise distributions since ESR currently supports only functions of one variable and so that we can evaluate the likelihood. Potential univariate problems in common benchmarks [14] are too simple or lack noise specifications. Problems 1 to 6 were created by hand by combining randomly sampled ESR solutions. Problems 7 to 12 were created by letting Operon with the same hyperparameters as in Table 1 with a maximum length of 50 over-fit to ten points of Gaussian random noise with zero mean and σ = 1 and then selecting and simplifying the best-fitting continuous function as ground truth. The values of x are evenly spaced in the given range to prevent bias from specific training or test set samples. The value range of x was chosen so that all values f (x) are on a similar scale. Small Gaussian noise with zero mean and σ = 10−7 is added to the target variable f (x) to focus on the optimization process rather than model selection. To verify the findings on synthetic data, we also run experiments on the onedimensional Nikuradse dataset [20], which contains measurements of friction of fluids in pipes [9]. Given the dense set of observations and its difficulty for modeling, this dataset has been used for testing SR implementations [9, 13] and allows us to verify whether the results hold in a more realistic setting. It contains 362 observations, which we randomly split into a training set of 230 and a test set of 132 observations [13]. Although the noise level of the Nikuradse dataset is not explicitly specified, it is expected to be much higher than the noise level of the synthetic problems. In our experiments, we assume a noise level of σ = 0.0367
6
L. Kammerer et al.
Table 2. Ground truth equations, their complexity and value ranges for the twelve synthetic problems. Gaussian noise N (0, 10−7 ) is added to f (x). Nr. Complexity Ground truth f (x) 1
11
2
13
3 4
13 11
x∈
1
[0.6, 2.5]
x x+x0.05 −x
2
−0.4x x − 0.5 + x
x
[0.1, 1.7]
−1 0.16 x xx
[0.6, 3]
x+3 x−x3.5
x
+ x1
13
0.1x +
6
17
7
41
x 0.1x + 0.10.17 − 0.05 x +0.3 x x 0.05 3 0.5x (−0.05x +(0.5x) )(0.1x+( 0.10.17 x +0.3 ) − x )
8
43
9
39
5
10 11
12
39 39
37
0.17 0.1x +0.25
[0.5, 2.5] x
[0.5, 2.5]
[0.2, 4]
− 0.06 0.15x+0.98 4.54 0.0061 |x − 1| + 0.348 + −0.02 + x 0.66 3.9 x 0.04x + 0.4 x − 1.2 + |x − 6.7|−2.1 x − 0.14 − 0.08 x 0.339 0.731.376
[0.08, 4.2] [0.5, 4]
1.661−x 1.095x 1.3760.36x (x−1.7) xx−1.7 −1.7 x−1.8 x −1.7
9.87 x −9.9x −x x0.723 x−0.72x (x(0.932x) ) x6.283 0.943x xx
0.47x−0.745 ) 0.636 (x + 0.318)(0.05x 0.745−0.47x 0.64 0.05
[0.5, 3]
[0.5, 1.5]
2 −x 1.4 0.05x −x
x−1.6 (x − 1.6) x0.049 + 0.27x0.073 − 0.028 0.104x−1.6 (2x−1.6) − 1.6 x0.036
−
[0.5, 3.2]
[0.2, 2]
determined as the root mean squared error (RMSE) of well-fitting models found with some exploratory GP runs. We use the same settings shown in Table 1, with a maximum complexity of 50.
3
Results
Figure 2 compares the best models found by ESR to the ground truth for each synthetic problem. Most ESR models capture the overall curvature, especially in Problems 1, 3 and 11. However, Problems 6, 7, and 12 show poor approximations. We report the root mean squared error (RMSE) in all our results which enables comparisons to the known irreducible noise level. The RMSE of the best ESR models ranges between 10−2 and 10−1 in all synthetic problems, which is much higher than the noise level of σ = 10−7 , despite visual similarity of the outputs. As population diversity is crucial for the success of GP and ESR does not explicitly enforce dissimilarity within its search space, we first analyze whether
Evaluation of Population Initialization Methods Best ESR Solutions fˆi (x) f (x), fˆi (x)
Problem 1
Ground Truth f (x) Problem 3
8
Problem 4 2
6
0.80 0.45
4 1
0.40
2
0.75 1
0.64
f (x), fˆi (x)
Problem 2
0.50
7
2
0.5
Problem 5
1.0
1.5
1
Problem 6
0.49
2
3
1
Problem 7
0.3
2
Problem 8 0.25
0.62
0.48 0.2
0.20
0.60 0.47 0.58 2
2
Problem 9 f (x), fˆi (x)
0.15
0.1 1
4
Problem 10
1.0
2
4
Problem 11
2 0.5
4
Problem 12
0
0.4 0.5
0.2
0
−2
0.0
−4
−0.5
0.0 1
2
x
3
0.5
1.0
x
1.5
1
2
x
3
1
2
x
Fig. 2. Ground truth and 100 most accurate ESR solutions. In most problems, ESR finds models that capture the overall curvature of the ground truth but still have a high RMSE of between 10−2 and 10−1 compared to the noise level of σ = 10−7 .
a population initialized with ESR provides sufficient diversity. To ensure that the subsequent evolutionary optimization is not compromised by a lack of diversity, we compare the diversity of an ESR-initialized population with randomly initialized populations using BTC, PTC-2 and the grow method. We use the hash-based tree distance [5] as a measure of syntactic dissimilarity between two expression trees. The hash-based tree distance is the ratio of distinct subtrees to the total number of subtrees in both expression trees, and also accounts for commutativity of multiplication and addition. Therefore, two trees that are equal or isomorphic regarding commutativity have a tree distance of zero, while two trees that do not share a single node have a tree distance of one. Figure 3 shows the mean and standard deviation of pairwise tree distances between the models in the populations of each run with error bars for each method. Since the ESRinitialized population depends on the problem, we show the distribution for each problem separately. For each random population initialization method, we show the mean across all 1000 runs as well as the average standard deviation of each of the 1000 populations. Since the random initial population is independent of the problem, we show it only once. The ESR-initialized populations show similar diversity across all problems. Their diversity is slightly below the diversity of populations initialized with BTC and PTC-2, but clearly within the standard deviation. The grow method shows slightly less diversity than the other methods with higher standard deviations,
8
L. Kammerer et al.
which is expected as it creates trees of any size and shape. Overall, the diversity of the ESR-initialized population is comparable to (although slightly lower than) the randomly initialized populations, suggesting that although the ESRinitialized population is more accurate, it is marginally less diverse which could hinder the evolutionary process. To quantify the accuracy advantage of the ESR-initialized population over the random initialization methods, we compare the distributions of test errors of models in the initial populations before any evolutionary optimization. Figure 4 shows the median and the centered 90% of test errors for each initialization method and problem. Invalid outputs that can occur in randomly generated models are treated as infinite. For the random initialization methods, we report the average of the median and the 5th and 95th percentiles over all 1000 runs. As expected, the ESR-initialized population achieves a median test RMSE one to two orders of magnitude lower than the random initialization methods, while still being clearly higher than the corresponding noise levels. The three random initialization methods yield very similar error distributions. To compare the error distributions of the final GP results of all four methods, we take the models from the Pareto front after 500 generations of all 1000 GP runs for each method and problem. Figure 5 shows the distribution of test error over complexity for each problem. The area between the dashed lines denotes the centered 90% of the test error of models at a specific complexity across all 1000 runs for a single method. The solid lines denote the median test error, and each dot indicates a specific complexity value. Given that ESR results are optimal with respect to accuracy within its search space, the dot at complexity 9 for GP with an ESR-initialized population shows the accuracy of the best ESR model.
Problem 1 Problem 2 Problem 3 Problem 4 Problem 5 Problem 6 Problem 7 ESR initialization Problem 8 Problem 9 Problem 10 Problem 11 Problem 12 Nikuradse BTC random initialization PTC-2 random initialization Grow random initialization 0.0
0.2
0.4
0.6
0.8
1.0
Tree Distance
Fig. 3. The mean and standard deviation of pairwise tree distances among individuals in the initial population for each method. Higher values denote more diversity. As ESR initialization depends on the problem, the mean and standard deviation are shown for each problem separately. For the three random initialization methods, the bars show the mean and the average standard deviation of the initial populations of all 1000 runs. All methods provide similar diversity, with just slightly different averages and the grow method with higher standard deviations.
Evaluation of Population Initialization Methods ESR initialization BTC random initialization
9
PTC-2 random initialization Grow random initialization
Problem 1 Problem 2 Problem 3 Problem 4 Problem 5 Problem 6 Problem 7 Problem 8 Problem 9 Problem 11 Problem 10 Problem 12 Nikuradse 10−2
10−1
100
101
102
RMSE
Fig. 4. The distribution of test errors of models in the initial population for all methods and problems. The dots show the median and the error bars denote the 5th and 95th percentile. The bars for the three random initialization methods show the average median and percentile values across all 1000 runs. As expected, we observe a much lower median test error in the ESR-initialized population than in the random initialization methods.
In all problems, this RMSE is between 10−2 and 10−1 . We omit the training error results as they are nearly identical to the shown test results. Despite the observed differences in the initial population with respect to diversity and accuracy, all methods perform nearly identically after the evolutionary process, with largely overlapping error distributions. Even the grow method with its slightly lower diversity and higher standard deviation performs equally well as the other methods. The ESR-initialized population also shows no clear improvement over GP with random population initialization and even performs worse in Problems 6 and 7. A visible improvement occurs only in Problem 1, where we reliably identify the ground truth with an ESR-initialized population. In most problems, such as Problem 11, all error distributions overlap completely. While the ESR-initialized population provides slightly better results after approximately the first twenty generations, the randomly initialized populations catch up quickly. Figures 6 and 7 show the test error distributions of models from the Pareto fronts of all runs after different generation counts for each initialization method for Problems 11 and 1, respectively. Figure 6 shows for Problem 11 that after generation one, GP with ESR initialization shows better results, but after ten generations, both error distributions largely overlap, with the ESRinitialized population showing only slightly better results. This difference shrinks further after 20 generations. After 200 generations, both methods show identical performance. This pattern is observed in all problems except Problem 1, which we discuss in more detail below. Many models created by the grow method ex-
10
L. Kammerer et al. ESR init., median BTC random init., median PTC-2 random init., median Grow random init., median Length of Ground Truth 10-2
Problem 1
RMSE
10-3
10
10-4
10
10-6
10-7 7
9 11 13 15 17 19
-2
Problem 5
RMSE
10-3
10
10-3 10-4
-5
-5
10
10
10-6
10-6
10-7
10-7
-2
7
9 11 13 15 17 19
Problem 9
10-3 10-5 10
Problem 6
10-5 10-7 7
10
-1
10
-3
10
-5
9 11 13 15 17 19
Problem 7
7 10-2
9 11 13 15 17 19
Problem 8
10-3
10-5 10-6
Problem 10
10-7 10-2
9 15 21 27 33 39 45
Problem 11
9 15 21 27 33 39 45 10
-1
Problem 12
10-2
10-4
10-3
10-5
-6
10-5
10-4
9 15 21 27 33 39 45 10-1
10-3
9 11 13 15 17 19
10-3
10-4
Problem 4
10-1
10-3
10-7 7
-2
10-4
10
RMSE
-5
10-6
10-7
Problem 3
10-1
-3
10-4
10-5
10
Problem 2
10-2
ESR init., 5th − 95th percentile BTC random init., 5th − 95th percentile PTC-2 random init., 5th − 95th percentile Grow random init., 5th − 95th percentile
10-6
10-4
9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
Complexity
Complexity
Complexity
Complexity
Fig. 5. The distribution of test error over complexity of the Pareto-optimal models from each of the 1000 runs per initialization method for the synthetic problems. The areas between the dashed lines denote centered 90% of test errors and the solid line the median at a specific complexity. All methods show similar performance with very similar distributions denoted by largely overlapping areas between the dashed lines. ESR initialization shows slightly better performance only in Problem 1 but performs even slightly worse in Problems 6 and 7.
ceed the complexity limits. Such models are not explicitly handled in Operon and appear on the Pareto front in the first few generations until dominated during the evolutionary process. We removed these models in early generations from our analysis as they would distort the results. We note that ESR initialization, which takes roughly two hours for all complexities per problem on an AMD EPYC 7713P 2.0 GHz 64-core CPU, is computationally much more expensive than the subsequent GP run, which took around only five seconds on the same CPU. These differences in computational effort outweigh any slight improvement in runtime of the GP run. The only clear improvement in accuracy and complexity by the ESR initialization is observed in Problem 1. We attribute this to the ESR solution 1 1 being very close to the ground truth f (x) = x+x0.05 fˆ(x) = x+x0.0974−x x −x , and therefore only requiring minor modifications by the evolutionary process. This is
Evaluation of Population Initialization Methods ESR init., median BTC random init., median PTC-2 random init., median Grow random init., median
11
ESR init., 5th − 95th percentile BTC random init., 5th − 95th percentile PTC-2 random init., 5th − 95th percentile Grow random init., 5th − 95th percentile Generation 20 Generation 200
Generation 1
Generation 10
9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
Complexity
Complexity
Complexity
Complexity
10-1
RMSE
10-2 10-3 10-4 10-5 10-6
Fig. 6. The distribution of test error over complexity of the Pareto-optimal models from each of the 1000 runs for each initialization method in the synthetic Problem 11 after four different numbers of generations. While the ESR-initialized population performs better initially, these differences shrink continuously after ten and 20 generations until they perform equally after 200 generations.
consistent with the observed search progress by GP shown in Figure 7, where GP with an ESR-initialized population identifies the ground truth after only around 20 generations in most runs. In contrast, GP with random initialization methods identifies the ground truth much less frequently even after 500 generations as shown in Figure 5. To ensure results are not affected by the specific algebraic representation of models selected by ESR, which uses symbolic simplification to prevent duplicates [3], we also run experiments with randomly sampled algebraic representations of the selected ESR models. These experiments show nearly identical results, with no significant improvement of the ESR-initialized population over random initialization in most problems. The results on the real-world Nikuradse dataset in Figure 8 show the same pattern as in the synthetic problems and verify that this behavior is also consistent in a setting containing more noise. Figure 8a plots the whole dataset as well as the curves of the 100 most accurate ESR models. It shows that the best ESR models provide only rough estimates with RMSE values clearly above the assumed noise level of σ = 0.03672. Figure 8b shows that all methods provide similar results after running GP for 500 generations with largely overlapping error distributions of the final Pareto fronts. Figure 9 also shows that, for this dataset, the initial advantage of an ESRinitialized population quickly diminishes after a few generations and does not lead to more accurate or less complex models. The same observation applies to all three random initialization methods, which perform equally well.
12
L. Kammerer et al.
ESR init., median BTC random init., median PTC-2 random init., median Grow random init., median Generation 1
10-1
RMSE
10
Generation 10
ESR init., 5th − 95th percentile BTC random init., 5th − 95th percentile PTC-2 random init., 5th − 95th percentile Grow random init., 5th − 95th percentile Generation 20 Generation 200
-2
10-3 10-4 10-5 10-6 7
9 11 13 15 17 19
7
Complexity
9 11 13 15 17 19
7
Complexity
9 11 13 15 17 19
7
Complexity
9 11 13 15 17 19
Complexity
Fig. 7. The distribution of test error over complexity of the Pareto-optimal models from each of the 1000 runs for each initialization method in the synthetic Problem 1 after four different numbers of generations. In contrast to all other problems, the ESR-initialized population shows a clear improvement over the randomly initialized population, which is due to the best ESR model being very close to the ground truth.
Training data Best ESR solutions fˆi (x)
ESR init. median/5th − 95th percentile BTC random init. median/5th − 95th percentile PTC-2 random init. median/5th − 95th percentile Grow random init. median/5th − 95th percentile
RMSE
y
2.0 1.5 1.0
9 × 10-2 8 × 10-2 7 × 10-2 6 × 10-2 5 × 10-2 4 × 10-2
0.5 0
1
2
3
9
15 21 27 33 39 45
x
Complexity
(a) The Nikuradse data and the 100 most accurate ESR models. The best ESR model achieves an RMSE of 0.077, around twice as high as the assumed noise level of σ = 0.03672.
(b) The distribution of test error over complexity of the Pareto-optimal models from each of the 1000 runs per initialization method [20, 9]. The areas between the dashed lines denote centered 90% of test errors and the solid line denotes the median at a specific complexity across all 1000 runs.
Fig. 8. The one-dimensional real-world Nikuradse dataset [20, 9] and the GP results for each initialization method, in which we observe a similar pattern as in the synthetic problems. While ESR-initialization provides reasonable fits in the initial population, as shown in Figure 8a, the GP yields similar results over all initialization methods, as denoted by the overlapping error distributions in Figure 8b.
Evaluation of Population Initialization Methods ESR init., median BTC random init., median PTC-2 random init., median Grow random init., median Generation 1
Generation 10
13
ESR init., 5th − 95th percentile BTC random init., 5th − 95th percentile PTC-2 random init., 5th − 95th percentile Grow random init., 5th − 95th percentile Generation 20 Generation 200
0.045
0.045
0.045
0.040
0.040
0.040
0.035
0.035
0.035
RMSE
0.20 0.15 0.10 0.05 9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
9 15 21 27 33 39 45
Complexity
Complexity
Complexity
Complexity
Fig. 9. The distribution of test error over complexity of the Pareto-optimal models from each of the 1000 runs for both methods in the Nikuradse dataset [20, 9] after four different numbers of generations. As in the synthetic problems, the advantage of the ESR-initialized population diminishes after only a few generations and does not lead to more accurate or less complex models.
4
Conclusion
Despite the intuitive assumption that optimizing the initial population for accuracy and complexity would improve the Pareto front of models in GP/SR, we find that an ESR-initialized population does not provide significant improvements over random initialization. Even with the best ESR population within the tested search space as GP’s initial population, GP with a randomly initialized population is able to recover any information provided by ESR on its own in the final Pareto front and quickly catches up with GP with an ESR-initialized population. Only in one problem, where the ground truth is very close to the best ESR model, do we observe an improvement. In all other problems, the initial advantage quickly diminishes within only a few generations. The same observation applies to the three random initialization methods, which all perform equally well. The grow method and its limited control over the complexity of the initial population led to models that exceed the given complexity limits and provided slightly less and more varying diversity. However, even these aspects do not impair the final results. This implies that the specific choice of initialization method barely has an effect on the final results in GP/SR. We recommend BTC or PTC-2 because they provide freely customizable complexity limits, are computationally inexpensive and lead to equivalent results in the subsequent evolutionary process as the grow method or ESR initialization. A limitation of our study is the small complexity in the optimized, ESRinitialized population. Although ESR provides optimal results within its supported complexity range, the error and complexity of the best ESR models are still far off from the noise level and the complexity of the ground truth. Since we observe a clear improvement in Problem 1 where the ground truth and the best ESR model are syntactically close, larger models in the initial population that
14
L. Kammerer et al.
are optimal in their covered search space might lead to more improvements over random initialization and increase the chance of finding a better or even optimal solution with GP. Even if such models in the larger-complexity initial population were not optimal, perhaps some ‘reasonable’ guesses based on domain-specific knowledge could be beneficial, as one would expect that minor modifications to these could achieve good solutions. We have only studied univariate problems. With higher-dimensional problems there may be differences between the random initialization methods. However, ESR cannot be used as it is limited to low dimensionality. Additionally, our results are based mostly on synthetic problems and only one real-world dataset. Other real-world problems might offer different possibilities for improvements, such as biasing the population towards specific functional forms based on prior knowledge about the underlying system. Finally, the complexity of the test functions also plays a role. We selected test functions with a range of complexities but potentially the datasets can be approximated well with short expressions which could make it easy for GP to evolve fitting expressions, even when the generating function is not recovered exactly. This work tests the effect of adapting only the initial population using ESR on the final Pareto front after many generations. Since we do observe beneficial effects of initialization, other approaches that maintain a stronger influence throughout the whole evolutionary process, such as periodic injection of elite individuals from ESR, might still lead to improvements over random initialization. Acknowledgement. L.K. and G.K. acknowledge support by the Austrian Federal Ministry for Economy, Energy, and Tourism, the Federal Ministry for Innovation, Mobility and Infrastructure, and the regional government of Upper Austria within the COMET project ProMetHeus (904919) supported by the Austrian Research Promotion Agency (FFG). D.J.B. acknowledges that support was provided by Schmidt Sciences, LLC. H.D. is supported by a Royal Society University Research Fellowship (grant no. 211046).
References 1. Ahmad, H., Helmuth, T.: A comparison of semantic-based initialization methods for genetic programming. In: Proceedings of the Genetic and Evolutionary Computation Conference Companion. pp. 1878–1881 (2018) 2. Bartlett, D.J., Desmond, H., Ferreira, P.G.: Exhaustive symbolic regression function sets (2022), https://doi.org/10.5281/zenodo.7339113 3. Bartlett, D.J., Desmond, H., Ferreira, P.G.: Exhaustive symbolic regression. IEEE Transactions on Evolutionary Computation 28(4), 950–964 (2023) 4. Burke, E., Gustafson, S., Kendall, G.: Ramped half-n-half initialisation bias in gp. In: Genetic and Evolutionary Computation Conference. pp. 1800–1801. Springer (2003) 5. Burlacu, B., Kammerer, L., Affenzeller, M., Kronberger, G.: Hash-based tree similarity and simplification in genetic programming for symbolic regression. In: International Conference on Computer Aided Systems Theory. pp. 361–369. Springer (2020)
Evaluation of Population Initialization Methods
15
6. Burlacu, B., Kronberger, G., Kommenda, M.: Operon C++ an efficient genetic programming framework for symbolic regression. In: Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion. pp. 1562–1570 (2020) 7. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T.: A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation 6(2), 182–197 (2002) 8. Deng, Y., Liu, Y., Zhou, D.: An improved genetic algorithm with initial population strategy for symmetric tsp. Mathematical problems in engineering 2015(1), 212794 (2015) 9. Guimerà, R., Reichardt, I., Aguilar-Mogas, A., Massucci, F.A., Miranda, M., Pallarès, J., Sales-Pardo, M.: A bayesian machine scientist to aid in the solution of challenging scientific problems. Science advances 6(5), eaav6971 (2020) 10. Hassanat, A.B., Prasath, V.S., Abbadi, M.A., Abu-Qdari, S.A., Faris, H.: An improved genetic algorithm with a new initialization mechanism based on regression techniques. Information 9(7), 167 (2018) 11. Koza, J.R.: Genetic Programming: On the Programming of Computers by Means of Natural Selection. MIT Press (1992) 12. Kronberger, G., Burlacu, B., Kommenda, M., Winkler, S.M., Affenzeller, M.: Symbolic regression. CRC Press (2024) 13. Kronberger, G., de Franca, F.O., Bartlett, D.J., Desmond, H., Ferreira, P.G.: Guiding multi-objective genetic programming with description length improves symbolic regression solutions. arXiv preprint arXiv:2605.22374 (2026) 14. La Cava, W., Orzechowski, P., Burlacu, B., de França, F.O., Virgolin, M., Jin, Y., Kommenda, M., Moore, J.: Contemporary symbolic regression methods and their relative performance. In: Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks. vol. 1 (2021) 15. Langdon, W.B., Nordin, J.: Seeding genetic programming populations. In: European Conference on Genetic Programming. pp. 304–315. Springer (2000) 16. Levenberg, K.: A method for the solution of certain non-linear problems in least squares. Quarterly of applied mathematics 2(2), 164–168 (1944) 17. Luke, S.: Two fast tree-creation algorithms for genetic programming. IEEE Transactions on Evolutionary Computation 4(3), 274–283 (2002) 18. Marquardt, D.W.: An algorithm for least-squares estimation of nonlinear parameters. Journal of the society for Industrial and Applied Mathematics 11(2), 431–441 (1963) 19. Mundhenk, T.N., Landajuela, M., Glatt, R., Santiago, C.P., Faissol, D.M., Petersen, B.K.: Symbolic regression via neural-guided genetic programming population seeding. In: Proceedings of the 35th International Conference on Neural Information Processing Systems. pp. 24912–24923 (2021) 20. Nikuradse, J.: Strömungsgesetze in rauhen Rohren. VDI-Forschungsheft 361, 1 (1933) 21. Paul, P.V., Ramalingam, A., Baskaran, R., Dhavachelvan, P., Vivekanandan, K., Subramanian, R.: A new population seeding technique for permutation-coded genetic algorithm: Service transfer approach. Journal of Computational Science 5(2), 277–297 (2014)