# Sharpness-Aware Minimization for Evolutionary Feature Construction in
Regression

Hengzhe Zhang, Qi Chen, Bing Xue, Wolfgang Banzhaf, Mengjie Zhang  H. Zhang,
Q. Chen, B. Xue and M. Zhang are with the School of Engineering and Computer
Science, Victoria University of Wellington, PO Box 600, Wellington 6140, New
Zealand (e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and mengjie.zhang@ecs.vuw.ac.nz).
W. Banzhaf is with the Department of Computer Science and Engineering,
Michigan State University, East Lansing, MI 48824, USA (e-mails:
banzhafw@msu.edu).

###### Abstract

In recent years, genetic programming (GP)-based evolutionary feature
construction has achieved significant success. However, a primary challenge
with evolutionary feature construction is its tendency to overfit the training
data, resulting in poor generalization on unseen data. In this research, we
draw inspiration from PAC-Bayesian theory and propose using sharpness-aware
minimization in function space to discover symbolic features that exhibit
robust performance within a smooth loss landscape in the semantic space. By
optimizing sharpness in conjunction with cross-validation loss, as well as
designing a sharpness reduction layer, the proposed method effectively
mitigates the overfitting problem of GP, especially when dealing with a
limited number of instances or in the presence of label noise. Experimental
results on 58 real-world regression datasets show that our approach
outperforms standard GP as well as six state-of-the-art complexity measurement
methods for GP in controlling overfitting. Furthermore, the ensemble version
of GP with sharpness-aware minimization demonstrates superior performance
compared to nine fine-tuned machine learning and symbolic regression
algorithms, including XGBoost and LightGBM.

###### Index Terms:

Automated Feature Engineering, Automated Machine Learning, Genetic
Programming, PAC-Bayesian, Generalization

##  1 Introduction

Automated feature construction is an important task in the domain of automated
machine learning (AutoML), which aims to improve the learning performance of a
machine learning algorithm on a given dataset $(X,Y)$ by constructing features
$\Phi(X)$ [1]. Automated feature construction methods include neural networks
[2], kernel methods [3], and genetic programming (GP). Among these, GP has
garnered significant attention in recent years due to its interpretable and
variable-length representation, as well as its gradient-free search mechanism,
making it well-suited for constructing complex and non-differentiable features
for non-differentiable base learners and loss functions [4].

However, the remarkable success of GP in feature construction is counteracted
by the challenge of overfitting [5], thereby constraining its broader
applicability [6, 7]. Current strategies for mitigating overfitting focus
primarily on either the structural risk minimization (SRM) framework or the
ensemble learning framework [8]. Examples of SRM-based approaches to combat
overfitting include the optimization of Vapnik-Chervonenkis (VC) dimension [6]
and Rademacher complexity [7]. On the other hand, ensemble learning methods
address overfitting by dynamically assembling an ensemble of top GP
individuals during the evolutionary process to construct features for decision
trees [9] and sigmoid functions [10].

However, a recent trend of viewing GP-based feature construction as an
efficient evolutionary deep learning technique [11] prompts us to rethink the
techniques for controlling overfitting in GP. In deep learning, the number of
parameters often exceeds the number of training instances, resulting in a high
VC dimension and Rademacher complexity for the neural network [12]. Despite
the VC dimension theory suggesting that large deep learning models can
memorize and thus overfit the training data [13], these deep models often
exhibit robust generalization performance. Therefore, relying on conventional
model complexity metrics to control overfitting in feature construction might
be too restrictive and insufficient, and alternative approaches should be
explored.

In the theoretical machine learning domain, there exists a theory known as
probably approximately correct (PAC)-Bayesian theory that has unveiled that
the expected generalization loss, denoted as $\mathbb{E}[\mathcal{L}(f)]$, for
a set of deep learning models $f$ with $n$ parameters following a normal
distribution $\mathcal{N}(\textbf{w},\bm{\sigma})^{n}$, is bounded as follows
with a probability of $1-\delta$ [14]:

$\begin{split}\mathbb{E}_{\bm{\nu}\sim\mathcal{N}(0,\sigma^{2})^{n}}\left[\mathcal{L}\left(f_{\mathbf{w}+\bm{\nu}}\right)\right]\leq&\underbrace{\mathbb{E}_{\bm{\nu}\sim\mathcal{N}(0,\sigma^{2})^{n}}\left[\widehat{\mathcal{L}}\left(f_{\mathbf{w}+\bm{\nu}}\right)\right]-\widehat{\mathcal{L}}\left(f_{\mathbf{w}}\right)}_{\text{Expected
Sharpness}}\\\
&+\underbrace{\widehat{\mathcal{L}}\left(f_{\mathbf{w}}\right)}_{\text{Training
Loss}}+\sqrt{\frac{1}{m}(\underbrace{\frac{\|\mathbf{w}\|_{2}^{2}}{2\sigma^{2}}}_{\mathrm{KL}}+\ln\frac{2m}{\delta})}\end{split}$
(1)

In this equation, $\widehat{\mathcal{L}}$ represents the empirical loss
function, $\bm{\nu}$ represents random perturbation parameters following a
normal distribution $\mathcal{N}(0,\sigma^{2})^{n}$, and $m$ denotes the
number of instances. Intuitively, the first term represents the change in loss
after a slight perturbation of parameters, the second term represents the
training loss, and the third term represents the deviation of parameters
$\mathbf{w}$ from a prior distribution $\mathcal{N}(0,\sigma^{2})^{n}$. Based
on this bound, the generalization loss of a machine learning model is
partially bounded by the change in loss when the model parameters are slightly
perturbed, which is known as expected sharpness.

To gain an intuitive understanding of the impact of sharpness minimization on
feature construction, i.e., minimizing the first term on the right-hand side
of Equation 1, Figure 1 illustrates an example of the training and test
$R^{2}$ scores using two different sets of constructed features on a dataset
named “Diabetes” [15]. One set consists of all possible combinations of third-
order addition features, such as $x_{1}+x_{2}+x_{3}$, as shown in Figure 2(a),
while the other comprises all possible third-order multiplication features,
like $x_{1}\cdot x_{2}\cdot x_{3}$, as depicted in Figure 2(b). In this case,
both the number of features and the size of features are identical, indicating
that model size cannot effectively distinguish the quality of these two groups
of features. Regarding training performance, multiplication features show a
better training $R^{2}$ score than addition features, as shown in Figure 1(a)
and Figure 1(c). However, on the test data, as illustrated in Figure 1(b) and
Figure 1(d), the additive features have a better test $R^{2}$ score than
multiplicative features. This suggests that additive features generalize much
better on the test data compared with multiplicative features. In this case,
model size cannot effectively control overfitting, but Figure 2(c) presents
promising results of using sharpness as an indicator to detect overfitting. In
Figure 2(c), we observe that the training loss of additive features is not
sensitive to perturbation, whereas the training loss of multiplicative
features is highly sensitive to noise/perturbation, implying that
multiplicative features correspond to minima of the training loss with high
sharpness. Therefore, although multiplicative features achieve higher accuracy
on the training data, they are not the preferred choice from the viewpoint of
sharpness-aware minimization (SAM) for achieving good generalization
performance.

(a) Additive Features

(b) Additive Features

(c) Multiplicative Features

(d) Multiplicative Features

Figure 1: Training and Test $R^{2}$ Scores for Different Constructed Features
on the “Diabetes” Dataset.

(a) Addition

(b) Multiplication

(c) Sharpness

Figure 2: Constructed Features and Their Sharpness Values.

Based on the above observations, this paper proposes a sharpness-aware
minimization framework for GP-based evolutionary feature construction
111Source Code:https://anonymous.4open.science/r/SAM-
GP/experiment/methods/SHARP_GP.py. Since GP optimizes the function space
rather than the parameter space, this paper proposes a sharpness estimation
method to gauge the sharpness of symbolic models. Then, both cross-validation
loss and the estimated sharpness are optimized within a multi-objective GP
framework to mitigate overfitting in GP. The main objectives of this paper are
summarized as follows:

  * •

Introducing a sharpness-aware minimization framework to the symbolic
optimization algorithm to address the issue of overfitting in GP-based
evolutionary feature construction methods.

  * •

Theoretically analyzing the proposed empirical sharpness estimation method for
symbolic models, such as GP, within the context of deep learning, considering
that GP is a data-efficient deep learning method.

  * •

Proposing a sharpness reduction layer in evolutionary feature construction to
explicitly minimize sharpness, thereby further enhancing generalization
performance on unseen data.

  * •

Proposing a bounded prediction strategy for robust predictions, assuming in-
distribution predictions, to further enhance generalization performance on
unseen data.

  * •

Proposing a least recently used caching technique to avoid the repeated
perturbation of identical features across different GP individuals, thereby
improving the efficiency of sharpness estimation.

  * •

Comparing the proposed algorithm with six other methods that optimize
different complexity measurements in GP to demonstrate the effectiveness of
sharpness-aware minimization, and with nine other machine learning algorithms
to demonstrate the effectiveness of GP-based learning methods.

The remaining sections of this paper are organized as follows. Section 2
provides an overview of related work on GP for feature construction and
overfitting control techniques. Section 3 presents the details of the proposed
method. Section 4 introduces the experimental settings. Section 5 shows
experimental results comparing the proposed method with other complexity
methods and other machine learning algorithms. Section 6 offers further
analysis regarding sharpness-aware minimization. Finally, in Section 7, we
provide conclusions drawn from the study and outline potential directions for
future research.

##  2 Related Work

###  2.1 Evolutionary Feature Construction

Evolutionary feature construction methods have undergone extensive research in
recent years. They can be classified into filter-based [1], wrapper-based [16,
17, 7], and embedded [18] methods based on the objective function used. In
filter-based feature construction methods, the constructed features are
evaluated without using any machine learning algorithm but using metrics, such
as impurity [1]. These approaches offer good time efficiency and generalize
well across different machine learning algorithms. However, the constructed
features may not achieve top-notch performance on a specific learning
algorithm. In contrast, wrapper-based feature construction methods construct
features tailored to a specific learning algorithm, thereby achieving good
performance on that the wrapped algorithm but at the expense of high
computational resources. Finally, embedded methods are as a compromise between
filter and wrapper methods, with a typical example being the symbolic
regression method [18], which embeds feature construction in the learning
process. In this work, we focus on the wrapper-based feature construction
method due to its superior learning performance [16, 17].

###  2.2 Controlling Overfitting in Genetic Programming

Overfitting is a critical concern within the GP community and has garnered
significant attention. Early efforts attempted to control overfitting by
managing the model size of GP models, but this approach has limited
effectiveness in many cases [19]. Thus, many methods have been developed to
improve the generalization performance of GP.

####  2.2.1 Statistical Learning-based Regularization

In recent years, model complexity metrics, such as the degree of the Chebyshev
polynomial approximation [20], VC dimension [6], and Rademacher complexity
[7], have been employed to enhance generalization performance of GP models.
These techniques have achieved notable success in controlling functional
complexity and improving generalization performance. In addition to model
complexity, other approaches from statistical learning theory, such as
residual analysis [21], robust loss function [22], and Bayesian model
selection [23], have also been explored for overfitting control. On the one
hand, traditional statistical learning theory often emphasizes model capacity,
implying that large/complex models that are capable of memorizing all training
data, such as generative pre-trained transformer (GPT), should not generalize
well. Specifically, theoretical work has proven that the lower bound of the VC
dimension for ReLU-based linear neural networks with $L$ layers and $W$
parameters is $\Omega(WL\log(W/L))$ [24]. However, on the other hand, much
state-of-the-art work has indicated that large deep learning models can often
generalize very well even with high VC dimension and Rademacher complexity
that might indicate overfitting of all training data [25].

####  2.2.2 Implicit Overfitting Control Methods

In parallel with the development of modern machine learning techniques, some
strategies like random sampling [26], semi-supervised learning [27], soft
target [28], modularization learning [29], concept drift detection [30] and
early stopping [31], have been employed to mitigate overfitting in GP. While
these techniques from the machine learning community are practically useful in
some cases, they lack a solid theoretical foundation from generalization
theory to explain the effectiveness of these overfitting control techniques.

####  2.2.3 Sharpness-based Regularization

Deep representation learning techniques have demonstrated that using a large
model for feature construction does not necessarily lead to serious
overfitting [12], indicating that the overfitting is more about the behavior
of a model, i.e., the semantics of the constructed features. Thus it is
necessary to further investigate overfitting control techniques for feature
construction. In the deep learning domain, PAC-Bayesian theory-based measures,
such as sharpness [32], are more effective in predicting the generalization
performance of deep learning models than classical complexity measures [33].
Thus, it becomes imperative to study overfitting control techniques based on
PAC-Bayesian theory and investigate their applicability to GP-based
evolutionary feature construction.

###  2.3 Sharpness Aware Minimization

The PAC-Bayesian theory has successfully provided practical bounds for the
generalization loss of deep neural networks [14], sparking a growing interest
in utilizing PAC-Bayesian bounds to enhance the generalization performance of
deep learning models. Based on the PAC-Bayesian bound presented in Equation 1,
the first term
$\mathbb{E}_{\bm{\nu}\sim\mathcal{N}(0,\sigma^{2})^{n}}\left[\widehat{\mathcal{L}}\left(f_{\mathbf{w}+\bm{\nu}}\right)\right]-\widehat{\mathcal{L}}\left(f_{\mathbf{w}}\right)$
has been widely employed to optimize to improve generalization performance,
which is known as sharpness-aware minimization (SAM) [32]. The underlying
intuition is to find a model situated at the flat minima of the loss
landscape, leading to a lower expected generalization loss. The SAM process
generally comprises two main steps:

  * •

Sharpness Maximization: Gradient ascent or random search is applied to perturb
the weights of a model $w$ within a certain radius $\rho$, aiming to find the
worst point of training loss in the neighborhood of the current model.

  * •

Loss Minimization: The gradient obtained from the perturbed model is utilized
to update the weights of the original model, guiding the model toward a flat
minimum.

To achieve sharpness-aware minimization effectively and efficiently in real-
world learning algorithms, several aspects have been taken into account and
studied:

  * •

Optimization Direction: Depending on whether the sharpness maximization
direction is computed based on the entire batch size $n$ or a subset of $m$
data samples, the approaches are categorized as $n$-sharpness and
$m$-sharpness. In practice, $m$-sharpness is more effective [34].
$m$-sharpness has a special case that focuses on maximizing sharpness around a
neighborhood $\rho$ independently for each data point, i.e., $m=1$. In this
case, optimizing $m$-sharpness can be equivalent to optimizing average
sharpness in PAC-Bayesian bounds when $\rho$ and the learning rate are
sufficiently small [35]. Moreover, during optimization, the gradient of the
original model may conflict with the gradient of the perturbed model.
Projecting the original gradient onto the perturbed gradient can yield better
results than optimizing either one alone [36].

  * •

Optimization Step Size: As weights may vary in scale, using adaptive
perturbation based on weight scale can enhance performance [37].

  * •

Efficient Perturbation: Even though estimating sharpness with gradient ascent
is computationally intensive, one-step gradient ascent can serve as an
approximation. To further reduce computational costs, sharpness can also be
estimated through random perturbation. This approach is effective if the
randomly perturbed weights are simultaneously optimized in parallel with the
original weights [38].

While sharpness-aware minimization has made significant advancements in the
domain of deep learning [32], it has not yet been explored in the context of
GP-based evolutionary feature construction. This is primarily because GP
optimizes the function space rather than the parameter space. Therefore, GP
lacks a clear definition of weights, and thus it is necessary to develop a
method for sharpness estimation that is suitable for GP, which is semantic
perturbation in this paper.

##  3 The Proposed Method

In this section, we introduce a framework for sharpness-aware minimization in
GP-based evolutionary feature construction. The section begins with an
overview of the sharpness-aware minimization framework for feature
construction. Subsequently, we present a practical method for estimating
sharpness in GP, supported by theoretical analysis. To further reduce the
sharpness of constructed features, a sharpness reduction layer is proposed.
Finally, we propose two strategies, namely bounded prediction and ensemble
learning, to enhance generalization performance further.

###  3.1 Algorithm Framework

The sharpness-aware feature construction algorithm proposed in this paper
adheres to the conventional framework of evolutionary feature construction
algorithms, as illustrated in Figure 3, and includes the following steps:

  * •

Initialization: Initially, $N$ individuals are generated randomly. Each
individual includes one GP tree initialized using the ramped half-and-half
method [39], representing a constructed feature.

  * •

Solution Evaluation: This phase consists of three steps.

    * –

Feature Construction: First, features $\Phi(X)$ are constructed based on the
GP trees $\Phi$.

    * –

Cross-Validation: Next, a linear model with L2 regularization is trained using
these constructed features $\Phi(X)$ to make predictions and calculate the
mean squared error on the training data, serving as the first objective
$O_{1}(\Phi)$. The leave-one-out cross-validation (LOOCV) scheme is employed
on the training set to mitigate overfitting. Specifically, for a
regularization strength $\alpha$, the weight $\bm{w}$ and leverage $\bm{h}$
are first calculated as shown in Equation 2:

$\begin{split}\bm{w}&=\left(\mathbf{X}^{\top}\mathbf{X}+\alpha\mathbf{I}\right)^{-1}\mathbf{X}^{\top}\mathbf{Y},\\\
\bm{h}&=\text{diag}\left(\mathbf{X}\left(\mathbf{X}^{\top}\mathbf{X}+\alpha\mathbf{I}\right)^{-1}\mathbf{X}^{\top}\right).\end{split}$
(2)

Then, the leave-one-out cross-validation errors on training set $\bm{e}$ are
computed based on the weight and leverage as shown in Equation 3:

$\bm{e}=\left(\frac{\mathbf{Y}-\mathbf{X}\bm{w}}{1-\bm{h}}\right)^{2}.$ (3)

    * –

Sharpness Estimation: Finally, the sharpness of the model is estimated,
serving as the second objective value $O_{2}(\Phi)$. Specifically, constructed
features are perturbed, and predictions are made based on perturbed features
to estimate sharpness. The method for estimating sharpness will be described
in Section 3.3.

  * •

Parent Selection: For a population $P$, parent individuals are selected based
on lexicase selection. Lexicase selection is a semantic-based selection
operator focusing on both improving accuracy and diversity, but does not
incorporate the sharpness objective. The sharpness objective is mainly
optimized through the survival selection operator after the offspring
generation stage. The core idea of lexicase selection is to iteratively
construct filters to eliminate individuals based on the median absolute
deviation of the loss $\epsilon_{k}$ for a randomly selected training instance
$k$, i.e., only keeping individuals $\Phi$ that satisfy
$\mathcal{L}_{k}(\Phi)\leq\min_{\Phi^{\prime}\in
P}\mathcal{L}_{k}(\Phi^{\prime})+\epsilon_{k}$, where $\min_{\Phi^{\prime}\in
P}\mathcal{L}_{k}(\Phi^{\prime})$ represents the best error on the instance.
This process is repeated until only one individual remains.

  * •

Offspring Generation: After parent individuals are selected, offspring are
generated by applying random subtree crossover and random subtree mutation
operators. Moreover, we also employ a tree addition operator and a tree
deletion operator to support adding and deleting features [16], allowing GP to
automatically determine the appropriate number of trees/features.
Specifically, the tree addition operator randomly generates a tree and adds it
to an individual, whereas the feature tree operator randomly selects a tree
and deletes it from an individual.

  * •

Survival Selection: Based on the two objectives, environmental selection
operators from NSGA-II [40] are used to select $|P|$ individuals from a
combination of $|P|$ parents and $|P|$ offspring. Briefly, the dominance
relationship and crowding distance between two selected individuals are
sequentially checked to determine the superior individual. Formally, an
individual $\Phi_{1}$ dominates another individual $\Phi_{2}$ if the following
two conditions are satisfied:

$\begin{split}\forall i\in\\{1,2\\},O_{i}(\Phi_{1})\leq O_{i}(\Phi_{2})\\\
\exists j\in\\{1,2\\},O_{j}(\Phi_{1})<O_{j}(\Phi_{2})\end{split}$ (4)

The crowding distance is defined as the average distance of a point to its two
nearest neighbors in the objective space. If a superior individual cannot be
identified using these criteria, an individual is randomly selected as the
parent.

  * •

Archive Maintenance: At this stage, the historically best-performing
individual based on the sum of cross-validation loss and the sharpness value,
i.e., $O_{1}(\Phi)+O_{2}(\Phi)$, is archived. Typically, only one individual
is stored, which is subsequently used for making predictions on unseen data.
However, more individuals could be stored in the archive to make an ensemble
prediction on unseen data, but an ensemble model may reduce interpretability.

Figure 3: Algorithm Framework

###  3.2 Sharpness Aware Minimization (SAM)

For any feature construction method, including neural networks and GP,
sharpness is defined as the increase in loss after a slight perturbation of
features $\Phi$. Assuming the perturbed features are denoted by
$\tilde{\Phi}$, the objective functions for SAM of a learning model $f$ are as
follows:

$\min_{\Phi}(\underbrace{\widehat{\mathcal{L}}(f(\Phi(X)),Y)}_{\text{Training
Loss}},\underbrace{\mathbb{E}_{\tilde{\Phi}}\;\widehat{\mathcal{L}}(f(\tilde{\Phi}(X)),Y)-\widehat{\mathcal{L}}(f(\Phi(X)),Y)}_{\text{Expected
Sharpness}})$ (5)

The first term is the objective traditionally optimized, while the second term
measures the sharpness of a set of features $\Phi$. For features $\Phi$
constructed by a neural network with a set of parameters $w$, sharpness is
often approximated using the one-step gradient ascent method [32]. Therefore,
sharpness is defined as
$\widehat{\mathcal{L}}\left(f(w+\rho\frac{\nabla_{w}L(x)}{\|\nabla_{w}L(x)\|_{2}},X),Y\right)-\widehat{\mathcal{L}}\left(f(w,X),Y\right)$,
where $\rho$ denotes the perturbation range, and
$\frac{\nabla_{w}L(x)}{\|\nabla_{w}L(x)\|_{2}}$ is the normalized gradient
used to determine the direction of perturbation. Thus, the empirical objective
functions become:

$\text{ minimize
}\left\\{\begin{array}[]{l}O_{1}(\Phi)=\widehat{\mathcal{L}}\left(f(w,X),Y\right)\\\
O_{2}(\Phi)=\widehat{\mathcal{L}}\left(f(w+\rho\frac{\nabla_{w}L(x)}{\|\nabla_{w}L(x)\|_{2}},X),Y\right)-O_{1}(\Phi)\end{array}\right.$
(6)

Here, $O_{1}(\Phi)$ is the empirical loss and $O_{2}(\Phi)$ is sharpness.
Regarding the perturbation term
$\rho\frac{\nabla_{w}L(x)}{|\nabla_{w}L(x)|_{2}}$, it can be calculated for
each instance or computed based on a batch of data, resulting in different
schemes known as 1-sharpness and n-sharpness. The idea of 1-sharpness or
n-sharpness can be applied to optimization in function space. However,
calculating the gradient $\rho\frac{\nabla_{w}L(x)}{\|\nabla_{w}L(x)\|_{2}}$
for GP-based feature construction is challenging due to the unclear definition
of learnable weights $w$ in function space. Therefore, we propose a method to
estimate sharpness in Section 3.3.

###  3.3 Empirical Sharpness Estimation

In this section, we introduce an empirical method to estimate sharpness in the
context of GP. The section begins by describing how to perturb the
semantics/outputs of GP trees to explore the neighborhood of GP trees in the
semantic space. Then, based on perturbed GP trees from Section 3.3.1, the
method for calculating the worst sharpness of a GP individual is introduced in
Section 3.3.2. This paper employs a 1-SAM strategy instead of an n-SAM
strategy to calculate the worst sharpness, and their differences in GP are
introduced in Section 3.3.3. The relationship between sharpness estimation in
GP and neural networks is illustrated in Section 3.3.4.

####  3.3.1 Semantic Perturbation

In contrast to NN, GP primarily focuses on optimizing the structure $\Phi$ in
function space rather than the weights $w$ in parameter space, making it
challenging to apply perturbations directly. To address this issue, we
introduce a method for perturbing the inputs of each non-terminal node
$\psi_{c}$, instead of perturbing the structure, i.e., perturbing the
semantics/outputs of each subtree of a GP tree, including terminal nodes.
Formally, suppose the semantics/outputs of a subtree in a GP tree $a$ are
$\psi_{a}(X)$ for input data $X$; then, the perturbed outputs
$\tilde{\psi}_{a}(X)$ are defined as:

$\tilde{\psi}_{a}(X)=\psi_{a}(X)+\psi_{a}(X)\cdot\mathcal{N}(0,\sigma^{2})$
(7)

where $\mathcal{N}(0,\sigma^{2})$ represents Gaussian noise generated with a
mean of zero and a standard deviation of $\sigma$. The Gaussian noise is
scaled by the semantics $\psi_{a}(X)$ to adjust the magnitude of the
perturbation to match that of the semantics. For a parent node
$\text{FunctionNode}_{c}$ in a GP tree that receives outputs from nodes $a$
and $b$ as inputs, the semantics of the GP tree $\psi_{c}$, based on a single
node $c$ and its perturbed subtree $\tilde{\psi}_{a}$ and $\tilde{\psi}_{b}$,
are defined by:

$\psi_{c}(X)=\text{FunctionNode}_{c}(\tilde{\psi}_{a}(X),\tilde{\psi}_{b}(X)).$
(8)

This equation illustrates perturbation on a single node. In a typical GP tree
with recursive subtrees, the perturbation process is recursively applied to
each subtree $\psi$ of a GP tree $\phi$, resulting in a perturbed feature
$\tilde{\phi}(X)$. All perturbed features constitute a perturbed individual
$\tilde{\Phi}(X)$.

####  3.3.2 Semantic Sharpness Estimation

The entire process for sharpness estimation is provided in Algorithm 1, which
consists of four stages:

  * •

Original Prediction (Lines 2-3): In the initial stage of sharpness estimation,
we obtain features constructed by GP trees, denoted as $\Phi(X)$, and the
corresponding predictions $\hat{Y}$ based on the trained linear model $LM$,
i.e., $\hat{Y}_{i}=LM(\Phi(X))$.

  * •

Cache-based Random Perturbation (Line 5): During each iteration of the
perturbation process, Gaussian noise $\mathcal{N}(0,\sigma^{2})$ is added to
the semantics of each subtree $\psi$ within each GP tree $\phi$ in the GP
individual $\Phi$, according to Equation 7 and Equation 8. The Gaussian noise
is controlled by the random seed $k\in K$. For this random perturbation
process, a least recently used cache is constructed for each random seed $k\in
K$. If a GP tree $\phi$ has been randomly perturbed in other GP individuals,
the perturbed feature $\tilde{\Phi}(X)$ is directly retrieved from the cache
to save computational resources.

  * •

Perturbed Prediction (Line 6): Once perturbed features $\tilde{\Phi}(X)$ are
generated, these features are fed into the trained linear model $LM$ to
produce perturbed predictions $\tilde{Y}$, i.e.,
$\tilde{Y}_{i}=LM(\tilde{\Phi}(X))$.

  * •

Sharpness Calculation (Lines 7-8): Finally, sharpness is calculated as the
difference between the mean squared errors of $\tilde{Y}$ and $\hat{Y}$,
specifically $(\tilde{Y}-Y)^{2}-(\hat{Y}-Y)^{2}$. Here, we use $1$-sharpness,
and thus, the worst mean squared error on one instance/fitness case is
recorded as $S_{i}$ separately.

To ensure the optimization of the worst-case sharpness [34] rather than the
average sharpness, random perturbations and sharpness calculations are
repeated $K$ times on each GP tree. The final sharpness value $S_{i}$ for each
instance $i$ is the maximum sharpness over $K$ iterations, i.e.,
$S_{i}=\max_{k=1}^{K}(\tilde{Y_{i,k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}$
across the $k$ rounds of perturbed predictions $\tilde{Y_{i,k}}$ for instance
$i$. For a GP individual, the sharpness value is the average sharpness value
across $N$ instances, i.e., $\sum_{i=1}^{N}S_{i}/N$. In summary, the objective
function of sharpness-aware minimization in GP is:

$\text{ minimize
}\left\\{\begin{array}[]{l}O_{1}(\Phi)=\sum_{i=1}^{N}(\hat{Y_{i}}-Y_{i})^{2}/N\\\
O_{2}(\Phi)=\sum_{i=1}^{N}S_{i}/N\end{array}\right.$ (9)

where $O_{1}(\Phi)$ represents the empirical mean squared error, and
$O_{2}(\Phi)$ is the estimated sharpness. In practice, $O_{1}(\Phi)$ is
replaced by the cross-validation error, as provided in Equation 3.

TABLE I: Strategies in Sharpness-Aware Minimization for Neural Networks and Evolutionary Feature Construction. All these strategies are integrated in this paper. Strategies | Neural Network | Evolutionary Feature Construction
---|---|---
Sharpness Estimation [32] | Gradient Ascent | Random Search
$m$-Sharpness [34] | Mini-batch Sharpness | Mini-batch Sharpness
Random Perturbation [38] | Random Perturbation on Weights | Random Perturbation on Semantics
Worst-case Perturbation [34] | Gradient Ascent | Maximum of multiple rounds Estimation
Adaptive Sharpness [37] | Perturbation Scaling According to Weights | Perturbation Scaling According to Semantics

####  3.3.3 1-SAM for GP

The sharpness estimation process for GP introduced in this paper is 1-SAM. An
example illustrating the difference between 1-SAM and n-SAM is shown in Figure
4. In the case of 1-SAM, the final sharpness is calculated based on the
maximum sharpness value $S_{i}$ within each instance $i\in\\{1,\dots,N\\}$
across $K$ rounds. Thus, the overall sharpness $S_{\text{1-SAM}}$ is defined
as:

$S_{\text{1-SAM}}=\frac{1}{N}\sum_{i=1,\dots,N}\max_{k=1,\dots,K}(\tilde{Y_{i}^{k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}.$
(10)

The mean of all the maximum sharpness values across all instances is computed
as the final sharpness, i.e., $\sum_{i=1}^{N}S_{i}$. In contrast, for n-SAM,
sharpness is calculated by perturbing the entire batch of data all at once,
rather than one by one. If there are $K$ rounds of perturbation, the largest
mean sharpness over the $K$ rounds is used as the final sharpness. In this
case, the overall sharpness $S_{\text{n-SAM}}$ is defined as:

$S_{\text{n-SAM}}=\max_{k=1,\dots,K}\frac{1}{N}\sum_{i=1,\dots,N}(\tilde{Y_{i}^{k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}.$
(11)

In this work, our goal is to optimize 1-sharpness instead of n-sharpness
because it has been shown to have superior performance in controlling
overfitting in deep learning [32].

Figure 4: Differences between 1-SAM and n-SAM. Values in the figure denote
sharpness for each training instance, with red blocks representing the final
selected sharpness values and blue blocks representing unused sharpness
values.

####  3.3.4 SAM for GP vs. NN

Table I provides a summary of the differences in strategies used in sharpness
estimation in NN and GP-based evolutionary feature construction. The symbolic
nature of GP and its non-differentiable properties lead to distinctions in the
sharpness estimation method. However, there is still a connection between
sharpness estimation in GP and NN. Figure 5 offers a dual perspective on the
sharpness estimation process proposed in this paper. From the viewpoint of a
GP tree, the sharpness estimation process perturbs the semantics of subtrees
within the tree. From the neural network perspective, the process perturbs the
inputs of each layer. In PAC-Bayesian theory [14], sharpness is defined as the
change in loss resulting from parameter perturbations. To clarify the
relationship between adding noise to weights and inputs at each layer,
theoretical analysis is presented in Section 3.4. This analysis demonstrates
that in the case of a neural network, adding noise to the inputs at each layer
is equivalent to perturbing the weights. Thus, in the context of GP, since it
is challenging to define weight perturbation in function space, this can be
achieved by perturbing the inputs of each subtree in GP, i.e., changing the
semantics/outputs of all the children of each subtree, rather than perturbing
the weights.

Figure 5: Two perspectives on sharpness estimation for GP: GP can be
considered an efficient deep learning technique [4, 41] when the GP tree is
rotated 90 degrees. Algorithm 1 Semantic Sharpness Estimation for GP Trees.

1:GP Tree $\Phi$, Inputs $X$, Target Outputs $Y$, Linear Model $LM$, Gaussian
Noise Standard Deviation $\sigma$, Number of Iterations $K$

2:Initialize sharpness $S\leftarrow 0$

3:$\Phi(X)\leftarrow$ Feature Construction ($X,\Phi$)

4:$\hat{Y}\leftarrow$ Prediction($LM,\Phi(X)$)

5:for $k=1,\ldots,K$ do

6: $\tilde{\Phi}(X)\leftarrow$ Feature Construction with Noise
($X,\Phi,\sigma$)

7: $\tilde{Y^{k}}\leftarrow$ Prediction($LM,\tilde{\Phi}(X)$)

8: for $i=1,\ldots,N$ do

9:
$S_{i}\leftarrow\max({S_{i},(\tilde{Y_{i}^{k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}})$

10:Sharpness Value $\sum_{i=1}^{N}S_{i}/N$

###  3.4 Theoretical Foundation of Sharpness Estimation

######  Theorem 1.

For a layer in a neural network represented as
$f(x)=g(\sum_{i=1}^{k}w_{i}x_{i})$ with an arbitrary activation function $g$,
adding adaptive noise $\mathcal{N}(0,w_{i}^{2})^{k}$ [14, 37] to weights $w$
is equivalent to adding noise of $\mathcal{N}(0,x_{i}^{2})^{k}$ to the inputs
of this layer, where $k$ denotes the dimensionality of the input data for this
layer, and $x_{i}$ represents the value for each dimension $i$ of input $x$.

Proof: Through mathematical transformation, the inputs of the activation
function, $\sum_{i=1}^{k}w_{i}x_{i}$, can be transformed as follows:

$\begin{split}\sum_{i=1}^{k}(w_{i}+\epsilon_{i})x_{i}&=\sum_{i=1}^{k}w_{i}(1+\epsilon_{i}/w_{i})x_{i}\\\
&=\sum_{i=1}^{k}w_{i}(x_{i}+x_{i}\epsilon_{i}/w_{i})\end{split}$ (12)

Given that $x_{i}$ and $w_{i}$ are constants when adding noise, let the noise
added to $x_{i}$ be denoted as $\theta_{i}$, i.e.,
$\theta_{i}=\epsilon_{i}w_{i}/x_{i}$. Then, $\theta_{i}$ follows the following
cumulative distribution:

$\begin{split}CDF(x_{i}\epsilon_{i}/w_{i})&=P(x_{i}\epsilon_{i}/w_{i}<k)\\\
&=P(\epsilon_{i}<kw_{i}/x_{i})\\\
&=\int_{-\infty}^{kw_{i}/x_{i}}\frac{1}{w_{i}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\epsilon_{i}}{w_{i}}\right)^{2}}d\epsilon_{i}\\\
&=\int_{-\infty}^{k}\frac{1}{w_{i}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{w_{i}\theta_{i}}{w_{i}x_{i}}\right)^{2}}\frac{w_{i}}{x_{i}}d\theta_{i}\\\
&=\int_{-\infty}^{k}\frac{1}{x_{i}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\theta_{i}}{x_{i}}\right)^{2}}d\theta_{i}\\\
\end{split}$ (13)

From the cumulative distribution, it is clear that
$\theta_{i}=x_{i}\epsilon_{i}/w_{i}$ follows a Gaussian distribution
$\mathcal{N}(0,x_{i}^{2})$. Therefore, adding noise to the inputs in each
dimension is equivalent to adding adaptive noise to the weight of each
dimension. Furthermore, since $k$ random noises are added to $k$ dimensions,
and $\epsilon_{i},\ldots,\epsilon_{k}$, are independently sampled, the
conclusion in one dimension can be extended to all dimensions. As a result,
Theorem 1 has been proven. An empirical verification is provided in the
supplementary material. Now, given that perturbing parameters is equivalent to
perturbing inputs, and considering that GP does not have a clear definition of
weights, it is reasonable to perturb the inputs of each layer in GP, i.e., the
semantics/outputs of the previous layer of each layer.

###  3.5 Sharpness Reduction Layer

The conventional GP operators, including mathematical operators such as
$\times$ and $\div$, are highly sensitive to even minor changes in input,
which can lead to altered computation results and result in non-zero
sharpness. To address this issue, we propose a sharpness reduction layer, the
key idea of which is to denoise corrupted results. During the training
process, the value of
$\psi^{*}_{c}(X)=\text{FunctionNode}_{c}(\psi_{a}(X),\psi_{b}(X))$ is stored
alongside node $c$. In the inference stage, a binary search is applied to the
pre-sorted $\psi_{c}(X)$ to search for the nearest value, $\psi_{c}^{+}(x)$,
that matches the subtree value $\psi_{c}(x)$ on the test instance $x$, as
outlined in Algorithm 2. Then, the subtree value on the test instance $x$,
$\psi_{c}(x)$, is replaced by the nearest value $\psi_{c}^{+}(x)$. For a GP
tree with multiple layers, the sharpness reduction layer is applied to each
layer to stabilize the semantics of constructed features. The sharpness
reduction layer avoids abrupt changes in constructed features in the
neighborhood of training samples, which could lead to sharp increases in
training loss. By applying this sharpness reduction layer, the final model
could be made more robust to unseen data and have improved generalization
performance.

Algorithm 2 Sharpness Reduction Layer

1:Training Data $\psi_{c}(X)$, Test Data $\psi_{c}(x)$

2:$\text{index}\leftarrow\text{Binary Search}(\psi_{c}(X),\psi_{c}(x))$

3:$\text{index}\leftarrow\text{clip}(\text{index},1,\text{len}(X)-1)$
$\triangleright$ Handle boundary values

4:$\psi_{c}(x_{\text{left}})\leftarrow\psi_{c}(X)[\text{index}-1]$

5:$\psi_{c}(x_{\text{right}})\leftarrow\psi_{c}(X)[\text{index}]$

6:if
$(\psi_{c}^{+}(x)-\psi_{c}(x_{\text{left}}))\leq(\psi_{c}(x_{\text{right}})-\psi_{c}^{+}(x))$
then

7: $\psi_{c}^{+}(x)\leftarrow\psi_{c}(x_{\text{left}})$

8:else

9: $\psi_{c}^{+}(x)\leftarrow\psi_{c}(x_{\text{right}})$

10:$\psi_{c}^{+}(x)$

###  3.6 Post-Processing Strategies

To further enhance the generalization performance of GP-based evolutionary
feature construction, we propose two strategies: bounded prediction and
ensemble learning.

####  3.6.1 Bounded Prediction

In machine learning, it is usually assumed that unseen data are generated from
the same distribution as the training data. Consequently, predictions should
not deviate significantly from the labels in the training data. To account for
this, during training, we record both the maximum target values
$\max_{i=1}^{N}Y_{i}$ and the minimum target values $\min_{i=1}^{N}Y_{i}$ from
the training dataset. Then, when making predictions on unseen data, we
constrain the prediction $\hat{Y_{i}}$ within the range defined by
$\min_{i=1}^{N}Y_{i}$ and $\max_{i=1}^{N}Y_{i}$. This constraint is commonly
found in many popular machine learning algorithms, such as decision trees and
k-nearest neighbors. Formally, the bounded prediction
$\hat{Y}_{\text{bounded}}$ is defined as:

$\hat{Y}_{\text{bounded}}=\begin{cases}\min_{i=1}^{N}Y_{i}&\text{if
}\hat{Y}<\min_{i=1}^{N}Y_{i}\\\ \max_{i=1}^{N}Y_{i}&\text{if
}\hat{Y}>\max_{i=1}^{N}Y_{i}\\\ \hat{Y}&\text{otherwise}\end{cases}$ (14)

By leveraging the assumption that the data are generated from the same
distribution, we anticipate that this bounded prediction step can enhance
generalization performance of the proposed method.

####  3.6.2 Ensemble Learning

Ensemble learning is an effective technique for reducing algorithm variance,
thereby improving generalization performance on unseen data. Given that GP is
a population-based optimization algorithm, it is natural to maintain an
archive $A$ of the top $|A|$ individuals based on the sum of the average
leave-one-out cross-validation loss and sharpness values during the
evolutionary process to form the ensemble. In other words, during the archive
maintenance stage in Section 3.1, not only the best individual is saved but
also the top $A$ individuals. Then, ensemble predictions are generated by
averaging the predictions from all models $\hat{Y}_{\Phi}$ within the
ensemble. Formally, the ensemble prediction $\hat{Y}_{\text{ensemble}}$ is
defined as follows:

$\hat{Y}_{\text{ensemble}}=\frac{1}{|A|}\sum_{\Phi\in A}\hat{Y}_{\Phi}$ (15)

It is worth noting that ensemble learning is an optional component for SAM-
based GP that can further enhance generalization performance. It is not
recommended when interpretability is the primary focus.

##  4 Experimental Settings

This section presents the experimental settings used to compare SAM and other
overfitting control techniques for GP-based feature construction, including
benchmark datasets, parameter settings, evaluation protocols, and baseline
methods.

###  4.1 Benchmark Datasets

In this paper, algorithms are validated using the Penn Machine Learning
Benchmark (PMLB) [42], which is a curated collection of datasets from OpenML
[43]. Our work focuses on real-world datasets. Thus, we exclude any datasets
generated by the Friedman, Feynman, or Strogatz functions. After applying
these exclusions, 58 datasets remain in the PMLB repository.

###  4.2 Parameter Settings

The parameter settings conform to the conventional settings of GP, as
summarized in Table II. The analytical quotient (AQ) [44] replaces the
division operator to prevent zero division errors. AQ is defined as
$AQ=\frac{a}{\sqrt{1+(b^{2})}}$, where $a$ and $b$ are input variables. Also,
the analytical log is used to replace the traditional log, which is defined as
$Log=log(1+\sqrt{x^{2}})$ for input $x$. The sum of the probabilities for
random mutation, random tree addition, and random tree deletion is set to 0.5
to explore the search space effectively [45]. Elitism is only applied in the
baseline algorithm when the environmental selection operator is not used, as
the environmental selection operator in NSGA-II inherently supports elitism.

TABLE II: Parameter settings for SAM-GP. Parameter | Value
---|---
Maximal Population Size | 200
Number of Generations | 100
Crossover and Mutation Rates | 0.9 and 0.1
Tree Addition Rate | 0.5
Tree Deletion Rate | 0.5
Initial Tree Depth | 0-3
Maximum Tree Depth | 10
Initial Number of Trees | 1
Maximum Number of Trees | 10
Elitism (Number of Individuals) | 1
Standard Deviation of Perturbation | 0.3
Iterations of Perturbation | 10
Functions |  +, -, *, AQ, Square, Sqrt, Abs, Log, Max, Min, Sin, Cos, Negative, Sigmoid

###  4.3 Performance Measure

To ensure the reliability of our conclusions, all experiments are conducted
over 30 independent runs, each initialized with a different random seed. In
each run, 100 samples are randomly selected from the dataset to be used as the
training data, while the remaining samples are employed as the test data. This
setup follows the pioneering work in studying the generalization of GP [46,
47] and is designed to simulate scenarios with limited data availability. For
datasets with fewer than 100 samples in either the training or test data, we
split these datasets at a ratio of 50:50 to avoid having too little training
and test data. For the experiments on label noise learning in Section 6.4, the
dataset is split into an 80:20 ratio for training and testing. Gaussian noise
$\mathcal{N}(0,1)$ is added to the labels of the training data. To eliminate
the influence of magnitude differences on the evaluation scores, we employ the
$R^{2}$ metric, which is defined as
$1-\frac{\sum_{i}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i}(y_{i}-\bar{y})^{2}}$, where
$\hat{y}_{i}$ denotes the prediction of the model for data instance $i$,
$y_{i}$ is the actual value for instance $i$, and $\bar{y}$ is the mean of the
actual values. Finally, a statistical comparison of the algorithms is
conducted using the Wilcoxon signed-rank test with a significance level of
0.05.

In this study, we investigate both single-model and ensemble feature
construction strategies. For single-model feature construction, the algorithm
is referred to as SAM-GP. As for ensemble feature construction, we set the
ensemble size to 100 and refer to the algorithm as SAM-EGP.

###  4.4 Baseline Methods

This paper evaluates the performance of SAM-GP against seven baseline GP
algorithms with various overfitting control methods:

  * •

Standard GP without Regularization (Standard GP): This is the standard GP
algorithm that uses leave-one-out cross-validation (LOOCV) loss as the sole
optimization objective.

  * •

Parsimony Pressure (PP) [48]: This approach uses the size of GP trees as a
measure of complexity to control overfitting.

  * •

Tikhonov Regularization (TK) [49]: This method employs zero-order Tikhonov,
i.e., the L2 norm of the outputs of the model, as a measure of complexity.
Zero-order Tikhonov is used because there is no significant difference between
using different orders of Tikhonov regularization for GP, and it is
computationally efficient [49].

  * •

Grand Complexity (GC) [49]: GC considers both model size and Tikhonov
complexity, using their dominance rank to determine overall complexity.

  * •

Rademacher Complexity (RC) [7]: This is a data-dependent complexity measure
that gauges the capacity of a machine learning model to fit a dataset.
Formally, for a set of models $\mathcal{F}$, with corresponding loss functions
denoted as $\mathcal{L}$, the Rademacher complexity is defined as follows:

$\operatorname{R}_{n}(\mathcal{L})=\mathbb{E}\left[\sup_{l\in\mathcal{L}}\frac{1}{n}\sum_{i=1}^{n}\zeta_{i}l\left(x_{i},y_{i}\right)\right],$
(16)

where $\zeta_{i}$ represents a randomly sampled Rademacher variable from
$\\{-1,1\\}$. For regression, the Rademacher complexity could be unbounded due
to the unbounded loss. Thus, for simplicity, we fit the machine learning model
based on constructed features using target $-\zeta_{i}y_{i}$ to obtain an
estimated Rademacher complexity. The loss function used in this paper is mean
squared loss, and thus when $\zeta_{i}=1$, $l\left(x_{i},y_{i}\right)$ will be
optimized toward $(2y_{i})^{2}$, and when $\zeta_{i}=-1$,
$l\left(x_{i},y_{i}\right)$ will be optimized toward 0.

  * •

Weighted MIC between Residuals and Variables (WCRV) [21]: WCRV aims to
minimize the high correlation between the chosen variables and residuals.
Here, the maximum information coefficient (MIC) is used as the metric to
measure the correlation between variables $x^{k}$ and residuals $R$.
Specifically, the WCRV of an individual $\Phi$ is defined as:

$\begin{split}\operatorname{WCRV}\left(\Phi\right)=&\sum_{\operatorname{MIC}_{x^{k},Y}>=mv}\operatorname{MIC}_{x^{k},Y}\times\operatorname{MIC}_{x^{k},R}\\\
&+\sum_{\operatorname{MIC}_{x^{k},Y}<mv}\left(1-\operatorname{MIC}_{x^{k},Y}\right).\end{split}$
(17)

For the selected variable with high MIC with target variable $Y$ that
surpasses the median value $mv$, i.e.,
$\sum_{\operatorname{MIC}_{x^{k},Y}>=mv}$, the MIC between variables $x^{k}$
and residuals $R$ is summed. For the unimportant variable, the negative MIC
between variables $x^{k}$ and target $Y$ is summed.

  * •

Correlation between the Input and Output Distances (IODC) [28]: IODC first
calculates pairwise Euclidean distances between inputs $D(x_{i},x_{j})$ and
pairwise Euclidean distances between outputs $D(y_{i},y_{j})$, where $x_{i}$
and $x_{j}$ represent two arbitrary input points, and $y_{i}$ and $y_{j}$
represent two arbitrary output points. This results in two pairwise distance
matrices for inputs and outputs, denoted as $I$ and $O$ respectively. Then, it
measures the Pearson correlation between these two distance matrices using the
formula:

$\operatorname{IODC}(\Phi)=\frac{\operatorname{Cov}(\mathrm{I},\mathrm{O})}{\sigma_{\mathrm{I}}\sigma_{\mathrm{O}}}$
(18)

where $Cov$ is the covariance matrix and $\sigma$ is the standard deviation.
In this equation, a high correlation indicates better smoothness and may be
less prone to overfitting.

All baseline methods follow the same framework as SAM-GP; the only difference
lies in replacing the sharpness objective with the respective metrics of these
methods. Bounded prediction is applied to all baseline methods to avoid
extreme outlier predictions. For these baseline methods, existing work often
selects the model with the best training error on the Pareto front as the
final model [7]. However, such a model can still be overly complex. Therefore,
this paper uses the minimum Manhattan distance (MMD) method [50] to determine
the knee point as the final model. First, two objectives are normalized based
on the best and worst objectives on the Pareto front, denoted as
$(\mathcal{O}_{1}(\Phi),\mathcal{O}_{2}(\Phi))$, since these complexity
measures do not have the same scale as mean square error. Then, MMD selects
the point with the best sum of the two normalized objectives, i.e.,
$\mathcal{O}_{1}(\Phi)+\mathcal{O}_{2}(\Phi)$, as the final model.

##  5 Experimental Results

This section presents the experimental results of GP using SAM and six other
overfitting control techniques, as well as the standard method for GP-based
feature construction on 58 datasets, and compares the test performance, the
training performance, the training time, and the model size.

###  5.1 Comparison on Test Performance

TABLE III: Statistical comparison of test $R^{2}$ scores when optimizing
various model complexity measures.

| PP | RC | GC | IODC | TK | WCRV | Standard GP
---|---|---|---|---|---|---|---
SAM | 32(+)/23($\sim$)/3(-) | 51(+)/7($\sim$)/0(-) | 37(+)/20($\sim$)/1(-) | 41(+)/16($\sim$)/1(-) | 50(+)/8($\sim$)/0(-) | 41(+)/13($\sim$)/4(-) | 36(+)/16($\sim$)/6(-)
PP | — | 36(+)/19($\sim$)/3(-) | 15(+)/36($\sim$)/7(-) | 23(+)/25($\sim$)/10(-) | 24(+)/33($\sim$)/1(-) | 20(+)/30($\sim$)/8(-) | 24(+)/26($\sim$)/8(-)
RC | — | — | 2(+)/25($\sim$)/31(-) | 4(+)/30($\sim$)/24(-) | 4(+)/32($\sim$)/22(-) | 8(+)/24($\sim$)/26(-) | 15(+)/10($\sim$)/33(-)
GC | — | — | — | 22(+)/28($\sim$)/8(-) | 25(+)/30($\sim$)/3(-) | 17(+)/37($\sim$)/4(-) | 24(+)/21($\sim$)/13(-)
IODC | — | — | — | — | 21(+)/23($\sim$)/14(-) | 19(+)/22($\sim$)/17(-) | 23(+)/14($\sim$)/21(-)
TK | — | — | — | — | — | 10(+)/27($\sim$)/21(-) | 10(+)/27($\sim$)/21(-)
WCRV | — | — | — | — | — | — | 18(+)/24($\sim$)/16(-)

####  5.1.1 General Analysis

The pairwise statistical comparison of test $R^{2}$ scores on 58 datasets is
presented in Table III, which compares six overfitting control methods with
standard GP that uses standard GP as the sole objective function. In general,
incorporating complexity measures can enhance generalization performance,
although the degree of improvement varies, and some methods also degrade
generalization on several datasets. Firstly, it is evident that SAM delivers
outstanding generalization performance, significantly outperforming standard
GP on 36 datasets and underperforming on 6 datasets out of 58, demonstrating
the advantage of using SAM for better generalization. As for other methods,
they outperform standard GP on at most 24 datasets, indicating that their
effectiveness is limited. Moreover, SAM outperforms other methods on more than
32 datasets, indicating that SAM is a better method for overfitting control in
general compared to other methods.

####  5.1.2 Sharpness vs. Model Size Control

One interesting phenomenon is that although PP is a poor overfitting control
method for GP-based symbolic regression [19], it remains a strong competitor
for controlling overfitting in GP-based evolutionary feature construction.
There are two potential reasons. First, previous research in symbolic
regression suggests that applying parsimony pressure can significantly reduce
model accuracy due to the small model size [7]. However, within the context of
evolutionary feature construction, several small and weak trees can still
produce accurate results if they complement each other [51]. Second, choosing
the “knee point” as the final model provides a good trade-off between
overfitting and underfitting, contrasting with traditional approaches that
select the model with the highest training accuracy as the final solution,
which has the risk of overfitting due to excessive model complexity. These
factors enable parsimony pressure to exhibit superior performance in
overfitting control. Nonetheless, SAM still outperforms it on 32 datasets,
while performing worse on only 3 datasets. Thus, we can conclude that
constructing simple features can reduce overfitting to some extent, but
considering the semantics of constructed features is a better way to control
overfitting.

####  5.1.3 Sharpness vs. Extreme Value Control

If a model outputs an extremely large value for a single training sample, even
if the mean squared error is small, the model may still be undesirable due to
non-smoothness. This led to the development of Tikhonov regularization for
overfitting control. However, the results show that SAM outperforms it on 50
datasets and is no worse on any dataset, indicating that Tikhonov
regularization is not an effective method for overfitting control. In fact,
the limited effectiveness of Tikhonov regularization in overfitting control
for GP was previously acknowledged [49]. The deficiency of Tikhonov
regularization led to the introduction of grand complexity, which combines
Tikhonov regularization with parsimony pressure. However, the results
demonstrate that although grand complexity outperforms Tikhonov regularization
in overfitting control, SAM still outperforms it on 37 datasets. Based on the
above analysis, we can conclude that controlling to avoid outputting extreme
values is not enough to achieve good generalization performance.

####  5.1.4 Evolutionary Plots

To gain a deeper understanding of the superior test performance of SAM, Figure
6 shows the evolutionary plots of test $R^{2}$ scores on four representative
datasets. The results show that using only LOOCV leads to overfitting of the
training data in early generations, resulting in decreased test $R^{2}$
scores, as observed in “OpenML_503” and ‘OpenML_522.” A similar situation can
also be seen when using IODC and PP, although the degree of overfitting has
been reduced. As for other methods, like GC and RC, they do not suffer from
the issue of severe overfitting, but their test $R^{2}$ scores are stuck at a
relatively low level. In contrast to these methods, SAM gradually improves
$R^{2}$ scores over the evolutionary process and finally surpasses all
baseline methods in the last generation. Combining these results, it is
evident that SAM is a promising method that exhibits good performance in
improving generalization on unseen data.

Figure 6: Evolutionary plots of test $R^{2}$ scores for different complexity
control methods.

###  5.2 Comparison on Training Performance

TABLE IV: Statistical comparison of training $R^{2}$ scores when optimizing
various model complexity measures.

| PP | RC | GC | IODC | TK | WCRV | Standard GP
---|---|---|---|---|---|---|---
SAM | 18(+)/19($\sim$)/21(-) | 47(+)/7($\sim$)/4(-) | 27(+)/24($\sim$)/7(-) | 25(+)/18($\sim$)/15(-) | 32(+)/13($\sim$)/13(-) | 27(+)/18($\sim$)/13(-) | 0(+)/4($\sim$)/54(-)
PP | — | 53(+)/5($\sim$)/0(-) | 36(+)/22($\sim$)/0(-) | 30(+)/22($\sim$)/6(-) | 33(+)/18($\sim$)/7(-) | 28(+)/22($\sim$)/8(-) | 0(+)/3($\sim$)/55(-)
RC | — | — | 1(+)/12($\sim$)/45(-) | 2(+)/13($\sim$)/43(-) | 1(+)/4($\sim$)/53(-) | 2(+)/19($\sim$)/37(-) | 0(+)/0($\sim$)/58(-)
GC | — | — | — | 11(+)/34($\sim$)/13(-) | 15(+)/26($\sim$)/17(-) | 13(+)/32($\sim$)/13(-) | 0(+)/0($\sim$)/58(-)
IODC | — | — | — | — | 18(+)/17($\sim$)/23(-) | 20(+)/17($\sim$)/21(-) | 0(+)/1($\sim$)/57(-)
TK | — | — | — | — | — | 25(+)/16($\sim$)/17(-) | 0(+)/5($\sim$)/53(-)
WCRV | — | — | — | — | — | — | 0(+)/4($\sim$)/54(-)

####  5.2.1 Overfitting of GP

To further confirm the presence of overfitting in GP, we present a pairwise
comparison of training $R^{2}$ scores in Table IV. The experimental results
indicate that optimizing using LOOCV as the sole objective leads to superior
training performance, outperforming SAM in 54 out of the 58 datasets. However,
when we analyze the training results in conjunction with the test results
presented in Table III, it becomes clear that relying solely on LOOCV results
in overfitting the training data. In contrast, SAM, which exhibited
significantly weaker training performance, achieved the best performance on
the test data. The evolutionary plots for training $R^{2}$ depicted in Figure
7 also support this observation. As illustrated, while optimizing LOOCV shows
improved $R^{2}$ scores throughout the evolution process, Figure 6 reveals
that these improved training $R^{2}$ scores correspond to decreased test
$R^{2}$ scores.

####  5.2.2 Alleviating Overfitting with Strong Regularization

When addressing severe overfitting, a common approach is to apply strong
regularization. However, as indicated by Figure 7, the use of strong
regularization significantly limits the training $R^{2}$ and adversely affects
generalization performance. For instance, optimizing using RC introduces a
strong regularization effect, resulting in notably lower training $R^{2}$
values compared to those achieved through optimizing LOOCV. However, this
method does not exhibit good generalization performance because of
underfitting. In other words, enhancing generalization requires careful
selection of the correct implicit regularization bias. Excessive
regularization, as demonstrated in these real-world problems, is
counterproductive to generalization performance.

####  5.2.3 Alleviating Overfitting by SAM

Figure 7 and Figure 6 provide a general view of the advantage of SAM over
other complexity measures. To further demonstrate that SAM can effectively
mitigate overfitting throughout the entire evolutionary process, Figure 8
illustrates the changes in training $R^{2}$ and test $R^{2}$ over the
evolutionary process. There is a growing trend in both training and test
$R^{2}$ scores as the generations increase. The Pearson correlations between
training $R^{2}$ and test $R^{2}$ are also presented in Figure 8 and Figure 9.
The results in Figure 8 show a correlation range of 0.85 to 0.99 between
training and test $R^{2}$ scores when using SAM, indicating a consistently
strong and positive correlation throughout the optimization process. In
comparison, the results in Figure 9 show a lower correlation between training
and test $R^{2}$ scores, indicating the presence of overfitting. All these
findings affirm that SAM is an effective strategy for optimizing GP in
scenarios where overfitting is a concern.

Figure 7: Evolutionary plots of training $R^{2}$ scores for various complexity
control methods. Figure 8: Evolutionary plots of training and test $R^{2}$
scores for SAM-GP. Figure 9: Evolutionary plots of training and test $R^{2}$
scores for GP without regularization.

###  5.3 Comparison of Training Time

In addition to predictive performance, running time is also an important
metric to consider when using a complexity measure for evolutionary feature
construction in real-world applications. In Figure 10, we present the
distribution of running times for various measures. The results indicate that
estimating sharpness is more time-consuming than complexity measures such as
model sizes and grand complexity. Optimizing model sizes takes a similar
amount of time compared to optimizing LOOCV, but optimizing SAM takes
approximately three times longer for training compared to optimizing LOOCV
alone or optimizing model sizes/grand complexity. This outcome was expected
because grand complexity calculations are based solely on the structure of GP
trees, whereas estimating sharpness requires calculating the change in the
semantics of GP trees under perturbations, which is more time-consuming.
Nonetheless, optimizing SAM is still more efficient than optimizing Rademacher
complexity. More importantly, considering the significant improvement in
generalization performance, the additional training time is deemed
justifiable. The overfitting problem cannot be solved by simply increasing the
number of generations in evolution, and it is still a fair comparison even
though the training time is larger than standard GP.

Figure 10: Distribution of training time across 58 datasets when optimizing
different complexity measures.

###  5.4 Comparison of Tree Size

In this section, we present comparisons of the tree size of GP-constructed
features. Although this aspect is not the primary focus of this paper, it is
relevant to interpretability and generalization performance. The historical
data related to each sharpness reduction layer is considered as one node.
First, as depicted in Figure 11, optimizing tree size or grand complexity
results in the smallest trees; however, optimizing sharpness does not achieve
comparable results because the tree size has not been explicitly optimized.
Compared to PP, SAM results in significantly larger trees in 41 datasets,
yields similar results in 10 datasets, and achieves better results in only 7
datasets. Although not comparable to PP, optimizing sharpness still moderately
controls the model size compared to standard GP, even though tree size is not
the direct target of optimization. Specifically, as shown in Table V, SAM
results in significantly smaller trees in 52 datasets compared to standard GP
and yields similar results in 6 datasets. This demonstrates that optimizing
sharpness can lead to more interpretable models compared to purely optimizing
training loss. Nonetheless, it is essential to emphasize again that smaller
models do not necessarily result in better generalization performance. For
example, optimizing grand complexity and parsimony pressure results in small
tree sizes but leads to poor generalization performance, as shown in Table
III. This further confirms existing findings in GP that a complexity measure
relying solely on tree size is insufficient [19], and semantics need to be
considered in the optimization to improve generalization performance.

Figure 11: Distribution of tree sizes across 58 datasets when optimizing
different model complexity measures. TABLE V: Statistical comparison of tree
sizes when optimizing different model complexity measures.

| PP | RC | GC | IODC | TK | WCRV | Standard GP
---|---|---|---|---|---|---|---
SAM | 7(+)/10($\sim$)/41(-) | 24(+)/19($\sim$)/15(-) | 10(+)/14($\sim$)/34(-) | 44(+)/6($\sim$)/8(-) | 46(+)/11($\sim$)/1(-) | 15(+)/17($\sim$)/26(-) | 52(+)/6($\sim$)/0(-)
PP | — | 54(+)/4($\sim$)/0(-) | 32(+)/20($\sim$)/6(-) | 57(+)/1($\sim$)/0(-) | 58(+)/0($\sim$)/0(-) | 44(+)/11($\sim$)/3(-) | 58(+)/0($\sim$)/0(-)
RC | — | — | 3(+)/5($\sim$)/50(-) | 32(+)/25($\sim$)/1(-) | 35(+)/21($\sim$)/2(-) | 13(+)/8($\sim$)/37(-) | 45(+)/12($\sim$)/1(-)
GC | — | — | — | 56(+)/2($\sim$)/0(-) | 58(+)/0($\sim$)/0(-) | 38(+)/17($\sim$)/3(-) | 58(+)/0($\sim$)/0(-)
IODC | — | — | — | — | 8(+)/42($\sim$)/8(-) | 5(+)/8($\sim$)/45(-) | 22(+)/32($\sim$)/4(-)
TK | — | — | — | — | — | 2(+)/9($\sim$)/47(-) | 29(+)/26($\sim$)/3(-)
WCRV | — | — | — | — | — | — | 54(+)/4($\sim$)/0(-)

###  5.5 Comparison with Other Machine Learning Algorithms

In addition to comparing different complexity measures, this paper
investigates the performance of various machine learning algorithms, following
the evaluation protocol of the state-of-the-art symbolic regression benchmark
[52], but using the same training set split method as employed in this paper.
The hyperparameters of all baseline algorithms are fine-tuned using cross-
validation on the training data before model fitting. The evaluated algorithms
encompass state-of-the-art machine learning methods, including XGBoost [53]
and LightGBM [54], as well as leading symbolic regression algorithms like
Operon [55]. Notably, due to the limited number of instances, hyperparameters
for these machine learning algorithms are tuned using a grid search instead of
successive halving grid search. The experimental results, including test
$R^{2}$ scores, model sizes, and training times, are presented in Figure 12.
The results on test $R^{2}$ scores reveal that SAM-GP not only outperforms
leading symbolic regression algorithms but also exceeds the performance of
well-known machine learning algorithms in terms of generalization performance
when data is limited. Figure 13 shows the pairwise comparison of different
algorithms using the Wilcoxon signed-rank test with Benjamini & Hochberg
correction. The results in Figure 13 demonstrate that the advantage of SAM-GP
over fine-tuned XGBoost is significant, indicating the superior generalization
performance of SAM-GP. Furthermore, as shown in Figure 12 and Figure 13, the
ensemble version of SAM-GP, i.e., SAM-EGP, significantly outperforms SAM-GP
and fine-tuned XGBoost, as well as other machine learning algorithms in terms
of test $R^{2}$ scores, further confirming the advantages of GP with SAM when
dealing with limited data. Regarding model size, SAM-GP is one order of
magnitude smaller than XGBoost, highlighting the superior interpretability of
GP-based evolutionary feature construction methods. As for training time, SAM-
GP takes longer than XGBoost but has a comparable training speed to Operon.
While training time remains an issue that needs addressing in the future,
given the superior predictive performance and interpretability of SAM-GP, it
is a tradeoff that is acceptable in many machine-learning scenarios.

Figure 12: Median $R^{2}$ scores, model sizes, and training time of different
algorithms on 58 regression problems. Figure 13: Pairwise statistical
comparisons of test $R^{2}$ scores on 58 regression problems.

##  6 Further Analysis

###  6.1 Analysis of Sharpness Reduction Layer

Figure 14 shows the impact of the sharpness reduction layer on the test
$R^{2}$ scores. For clarity, if improvements or decreases in $R^{2}$ are less
than 0.01, they are simply denoted as neutral results, and the analysis is
more focused on those datasets where $R^{2}$ could improve or decrease by more
than 0.1. The results indicate that the sharpness reduction layer can improve
$R^{2}$ by more than 0.01 on 14 datasets compared to not using the sharpness
reduction layer, and it leads to an $R^{2}$ decrease of more than 0.01 on only
6 datasets. These results demonstrate that the sharpness reduction layer is an
effective way to improve generalization performance.

Figure 14: Improvements and decreases of more than 0.01 in $R^{2}$ scores
brought by using the sharpness reduction layer, denoted by red and blue dots,
respectively. The green dot represents neutral results with improvements and
decreases of less than 0.01.

###  6.2 Comparisons on Caching Technique

The experimental results investigating the effectiveness of the proposed
caching technique for accelerating the learning process are presented in
Figure 15. When using the caching technique, the median training time is
1036s. In contrast, without using the caching technique, the median training
time is 1893s. In other words, after implementing the caching technique, the
median training time is reduced by 45%. These experimental results show that
the proposed caching technique can significantly reduce the time required for
sharpness estimation.

Figure 15: Distribution of training time (seconds) with and without the use of
the caching technique.

###  6.3 Analysis of Bounded Prediction

In this paper, we introduce a bounded prediction mechanism to further reduce
sharpness and improve generalization performance. As depicted in Figure 16(b),
employing this mechanism along with SAM significantly improves generalization
performance on 7 datasets without reducing performance on any dataset. This
suggests the effectiveness of bounded predictions in GP-based evolutionary
feature construction. Figure 16(a) presents a statistical test comparison on
training performance, revealing that the proposed bounded prediction strategy
does not significantly impact training performance, which aligns with the
intuition that bounded prediction mainly restricts unreasonable extrapolation
on unseen data. To further understand the impact of bounded predictions, a
case study on “OpenML 537” is depicted in Figure 17. This figure shows that
the model underestimates and overestimates predictions near the left and right
boundaries of the prediction interval, respectively. The clipping strategy
makes the model provide conservative predictions, thus achieving better
prediction performance. The effectiveness of this clipping strategy in
regression echoes a recent success in controlling overfitting by clipping
logits in classification models [56]. In the future, this issue could be more
effectively addressed by finding a method to measure the KL divergence of GP
models to control the KL divergence term in the PAC-Bayesian bound, i.e.,
$\frac{\|\mathbf{w}\|_{2}^{2}}{2\sigma^{2}}$, to guide the model in avoiding
predictions that deviate excessively from prior knowledge and further enhance
generalization performance.

(a) Training $R^{2}$

(b) Test $R^{2}$

Figure 16: Statistical comparison of training and test $R^{2}$ scores before
and after using bounded predictions. Figure 17: Effect of bounded predictions
on “OpenML 537” test data.

###  6.4 Performance in the Presence of Label Noise

In addition to limited samples, another scenario in evolutionary feature
construction that poses a risk of overfitting is label noise. In the domain of
deep learning, sharpness-aware minimization has demonstrated effectiveness in
addressing label noise [32]. In this section, we investigate the effectiveness
of the proposed SAM-GP in label noise scenarios. Due to limited computational
resources, the experiments are conducted on 36 datasets with a maximum of 2000
instances, and only standard GP and the top-performing method, i.e., PP, are
compared. The experimental results are shown in Figure 18. The results
demonstrate that SAM improves generalization performance compared to standard
GP on 31 datasets and outperforms PP on 21 datasets, without showing worse
performance on any dataset. These results indicate that using sharpness-aware
minimization in evolutionary feature construction is beneficial when labels
contain noise, consistent with the findings in the deep learning domain [32].

Figure 18: Statistical comparison of test $R^{2}$ scores on 36 datasets with
label noise.

###  6.5 Visualization of Constructed Features

In this paper, we have demonstrated the superior performance of SAM-GP. In
this section, we investigate the features constructed by SAM-GP and standard
GP to gain a deeper understanding of the differences between controlling
overfitting by optimizing sharpness or not. An example of features optimized
by standard GP and SAM-GP on “OpenML 522” is illustrated in Figure 19(a).
Experimental results reveal that the features constructed using standard GP
are complex, including nested functions like $sin(sin())$. In contrast, the
model obtained through SAM is less functionally complex. In terms of training
performance, standard GP can lead to a training performance of 0.22, while the
training $R^{2}$ of the model obtained through SAM is 0.16. However, regarding
test performance, SAM achieves a test $R^{2}$ score of 0.18, whereas the model
obtained by standard GP only has a test $R^{2}$ score of 0.0, indicating
overfitting in the model learned through standard GP and good generalization
performance of the model learned by SAM.

(a) No Regularization

(b) Sharpness-Aware Minimization

Figure 19: Examples of constructed features based on (a). no regularization
(b). sharpness-aware minimization.

##  7 Conclusions

To enhance the generalization ability of evolutionary feature construction,
this paper proposes a sharpness-aware minimization technique for symbolic
models to control the sharpness of GP-constructed features and generate
features that generalize well to unseen data. The experimental results on 58
datasets demonstrate the effectiveness of the proposed method in improving the
generalization performance of constructed features compared to solely
optimizing cross-validation loss or combining well-known complexity measures
in statistical machine learning. This improvement surpasses that achieved by
optimizing six other complexity measures. When comparing the optimization of
SAM with other metrics, we find that improving generalization is more about
achieving a flat minimum in the semantic space rather than optimizing model
size. Additionally, the experimental results indicate that SAM-GP
significantly outperforms fine-tuned state-of-the-art machine learning
algorithms on limited samples, further demonstrating the effectiveness of GP-
based learning algorithms.

In this work, SAM is not the best option for all datasets. Therefore, in the
future, it is worthwhile to determine the most suitable situations for using
SAM based on the characteristics of each dataset, such as employing the
instance space analysis technique [57]. Another potential future direction is
to validate the proposed technique with other AutoML approaches. The proposed
method is not limited to evolutionary-based methods but is also applicable to
Monte Carlo tree search-based AutoML methods [58] and Bayesian optimization-
based methods [59]. Finally, it is also worth exploring the applicability of
the proposed method to various symbolic regression algorithms, including
reinforcement learning-based [60] and Transformer-based [61] symbolic
regression techniques.

## References

  * Neshatian et al. [2012] K. Neshatian, M. Zhang, and P. Andreae, “A filter approach to multiple feature construction for symbolic learning classifiers using genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 16, no. 5, pp. 645–661, 2012.
  * LeCun et al. [2015] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _Nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * Liu et al. [2021] F. Liu, X. Huang, Y. Chen, and J. A. Suykens, “Random features for kernel approximation: A survey on algorithms, theory, and beyond,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 10, pp. 7128–7148, 2021.
  * La Cava et al. [2018] W. La Cava, T. R. Singh, J. Taggart, S. Suri, and J. H. Moore, “Learning concise representations for regression by evolving networks of trees,” in _International Conference on Learning Representations_ , 2018.
  * Agapitos et al. [2019] A. Agapitos, R. Loughran, M. Nicolau, S. Lucas, M. O’Neill, and A. Brabazon, “A survey of statistical machine learning elements in genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 23, no. 6, pp. 1029–1048, 2019.
  * Chen et al. [2018] Q. Chen, M. Zhang, and B. Xue, “Structural risk minimization-driven genetic programming for enhancing generalization in symbolic regression,” _IEEE Transactions on Evolutionary Computation_ , vol. 23, no. 4, pp. 703–717, 2018\.
  * Chen et al. [2022] Q. Chen, B. Xue, and M. Zhang, “Rademacher complexity for enhancing the generalization of genetic programming for symbolic regression,” _IEEE Transactions on Cybernetics_ , vol. 52, no. 4, pp. 2382–2395, 2022.
  * Nag and Pal [2015] K. Nag and N. R. Pal, “A multiobjective genetic programming-based ensemble for simultaneous feature selection and classification,” _IEEE Transactions on Cybernetics_ , vol. 46, no. 2, pp. 499–510, 2015.
  * Zhang et al. [2023a] H. Zhang, A. Zhou, Q. Chen, B. Xue, and M. Zhang, “Sr-forest: A genetic programming based heterogeneous ensemble learning method,” _IEEE Transactions on Evolutionary Computation_ , 2023.
  * Wei et al. [2020] T. Wei, W.-L. Liu, J. Zhong, and Y.-J. Gong, “Multiclass classification on high dimension and low sample size data using genetic programming,” _IEEE Transactions on Emerging Topics in Computing_ , vol. 10, no. 2, pp. 704–718, 2020.
  * Bi et al. [2022] Y. Bi, B. Xue, and M. Zhang, “Genetic programming-based evolutionary deep learning for data-efficient image classification,” _IEEE Transactions on Evolutionary Computation_ , 2022.
  * Zhang et al. [2021a] C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals, “Understanding deep learning (still) requires rethinking generalization,” _Communications of the ACM_ , vol. 64, no. 3, pp. 107–115, 2021.
  * Scarselli and Tsoi [1998] F. Scarselli and A. C. Tsoi, “Universal approximation using feedforward neural networks: A survey of some existing methods, and some new results,” _Neural Networks_ , vol. 11, no. 1, pp. 15–37, 1998.
  * Neyshabur et al. [2017] B. Neyshabur, S. Bhojanapalli, D. McAllester, and N. Srebro, “Exploring generalization in deep learning,” _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * Efron et al. [2004] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” _Annals of Statistics_ , pp. 407–451, 2004.
  * La Cava and Moore [2020] W. La Cava and J. H. Moore, “Learning feature spaces for regression with genetic programming,” _Genetic Programming and Evolvable Machines_ , vol. 21, pp. 433–467, 2020.
  * Zhang et al. [2022] H. Zhang, A. Zhou, and H. Zhang, “An evolutionary forest for regression,” _IEEE Transactions on Evolutionary Computation_ , vol. 26, no. 4, pp. 735–749, 2022.
  * Chen et al. [2017] Q. Chen, M. Zhang, and B. Xue, “Feature selection to improve generalization of genetic programming for high-dimensional symbolic regression,” _IEEE Transactions on Evolutionary Computation_ , vol. 21, no. 5, pp. 792–806, 2017\.
  * Vanneschi et al. [2010] L. Vanneschi, M. Castelli, and S. Silva, “Measuring bloat, overfitting and functional complexity in genetic programming,” in _Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation_ , 2010, pp. 877–884.
  * Vladislavleva et al. [2008] E. J. Vladislavleva, G. F. Smits, and D. Den Hertog, “Order of nonlinearity as a complexity measure for models generated by symbolic regression via pareto genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 13, no. 2, pp. 333–349, 2008.
  * Chen et al. [2020] Q. Chen, B. Xue, and M. Zhang, “Improving symbolic regression based on correlation between residuals and variables,” in _Proceedings of the 2020 Genetic and Evolutionary Computation Conference_ , 2020, pp. 922–930.
  * Raymond et al. [2023] C. Raymond, Q. Chen, and B. Xue, “Learning symbolic model-agnostic loss functions via meta-learning,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2023.
  * Bomarito et al. [2022] G. F. Bomarito, P. E. Leser, N. Strauss, K. M. Garbrecht, and J. D. Hochhalter, “Bayesian model selection for reducing bloat and overfitting in genetic programming for symbolic regression,” in _Proceedings of the Genetic and Evolutionary Computation Conference Companion_ , 2022, pp. 526–529.
  * Bartlett et al. [2019] P. L. Bartlett, N. Harvey, C. Liaw, and A. Mehrabian, “Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks,” _The Journal of Machine Learning Research_ , vol. 20, no. 1, pp. 2285–2301, 2019.
  * Jiang et al. [2019] Y. Jiang, B. Neyshabur, H. Mobahi, D. Krishnan, and S. Bengio, “Fantastic generalization measures and where to find them,” in _International Conference on Learning Representations_ , 2019.
  * Gonçalves and Silva [2013] I. Gonçalves and S. Silva, “Balancing learning and overfitting in genetic programming with interleaved sampling of training data,” in _Genetic Programming: 16th European Conference, EuroGP 2013, Vienna, Austria, April 3-5, 2013. Proceedings 16_. Springer, 2013, pp. 73–84.
  * Silva et al. [2018] S. Silva, L. Vanneschi, A. I. Cabral, and M. J. Vasconcelos, “A semi-supervised genetic programming method for dealing with noisy labels and hidden overfitting,” _Swarm and Evolutionary Computation_ , vol. 39, pp. 323–338, 2018.
  * Vanneschi and Castelli [2021] L. Vanneschi and M. Castelli, “Soft target and functional complexity reduction: A hybrid regularization method for genetic programming,” _Expert Systems with Applications_ , vol. 177, p. 114929, 2021.
  * Zhang et al. [2023b] H. Zhang, Q. Chen, B. Xue, W. Banzhaf, and M. Zhang, “Modular multi-tree genetic programming for evolutionary feature construction for regression,” _IEEE Transactions on Evolutionary Computation_ , 2023.
  * Celik and Vanschoren [2021] B. Celik and J. Vanschoren, “Adaptation strategies for automated machine learning on evolving data,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 43, no. 9, pp. 3067–3078, 2021.
  * Tuite et al. [2011] C. Tuite, A. Agapitos, M. O’Neill, and A. Brabazon, “Early stopping criteria to counteract overfitting in genetic programming,” in _Proceedings of the 13th Annual Conference Companion on Genetic and Evolutionary Computation_ , 2011, pp. 203–204.
  * Foret et al. [2020] P. Foret, A. Kleiner, H. Mobahi, and B. Neyshabur, “Sharpness-aware minimization for efficiently improving generalization,” in _International Conference on Learning Representations_ , 2020.
  * Owen et al. [2022] C. A. Owen, G. Dick, and P. A. Whigham, “Standardization and data augmentation in genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 26, no. 6, pp. 1596–1608, 2022.
  * Andriushchenko and Flammarion [2022] M. Andriushchenko and N. Flammarion, “Towards understanding sharpness-aware minimization,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 639–668.
  * Wen et al. [2022] K. Wen, T. Ma, and Z. Li, “How sharpness-aware minimization minimizes sharpness?” in _The Eleventh International Conference on Learning Representations_ , 2022.
  * Zhuang et al. [2021] J. Zhuang, B. Gong, L. Yuan, Y. Cui, H. Adam, N. C. Dvornek, J. s Duncan, T. Liu _et al._ , “Surrogate gap minimization improves sharpness-aware training,” in _International Conference on Learning Representations_ , 2021\.
  * Kwon et al. [2021] J. Kwon, J. Kim, H. Park, and I. K. Choi, “Asam: Adaptive sharpness-aware minimization for scale-invariant learning of deep neural networks,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 5905–5914.
  * Li et al. [2022] T. Li, W. Yan, Z. Lei, Y. Wu, K. Fang, M. Yang, and X. Huang, “Efficient generalization improvement guided by random weight perturbation,” _arXiv preprint arXiv:2211.11489_ , 2022.
  * Banzhaf et al. [1998] W. Banzhaf, P. Nordin, R. E. Keller, and F. D. Francone, _Genetic programming: an introduction: on the automatic evolution of computer programs and its applications_. Morgan Kaufmann Publishers Inc., 1998.
  * Deb et al. [2002] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: Nsga-ii,” _IEEE Transactions on Evolutionary Computation_ , vol. 6, no. 2, pp. 182–197, 2002.
  * Gaier and Ha [2019] A. Gaier and D. Ha, “Weight agnostic neural networks,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * Olson et al. [2017] R. S. Olson, W. La Cava, P. Orzechowski, R. J. Urbanowicz, and J. H. Moore, “PMLB: a large benchmark suite for machine learning evaluation and comparison,” _BioData mining_ , vol. 10, no. 1, pp. 1–13, 2017.
  * Vanschoren et al. [2014] J. Vanschoren, J. N. Van Rijn, B. Bischl, and L. Torgo, “Openml: networked science in machine learning,” _ACM SIGKDD Explorations Newsletter_ , vol. 15, no. 2, pp. 49–60, 2014.
  * Ni et al. [2012] J. Ni, R. H. Drieberg, and P. I. Rockett, “The use of an analytic quotient operator in genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 17, no. 1, pp. 146–152, 2012.
  * Muñoz et al. [2019] L. Muñoz, L. Trujillo, S. Silva, M. Castelli, and L. Vanneschi, “Evolving multidimensional transformations for symbolic regression with m3gp,” _Memetic Computing_ , vol. 11, pp. 111–126, 2019.
  * Ni and Rockett [2015] J. Ni and P. Rockett, “Training genetic programming classifiers by vicinal-risk minimization,” _Genetic Programming and Evolvable Machines_ , vol. 16, pp. 3–25, 2015.
  * Nicolau and Agapitos [2021] M. Nicolau and A. Agapitos, “Choosing function sets with better generalisation performance for symbolic regression models,” _Genetic programming and evolvable machines_ , vol. 22, no. 1, pp. 73–100, 2021.
  * Zhang and Mühlenbein [1995] B.-T. Zhang and H. Mühlenbein, “Balancing accuracy and parsimony in genetic programming,” _Evolutionary Computation_ , vol. 3, no. 1, pp. 17–38, 1995.
  * Ni and Rockett [2014] J. Ni and P. Rockett, “Tikhonov regularization as a complexity measure in multiobjective genetic programming,” _IEEE Transactions on Evolutionary Computation_ , vol. 19, no. 2, pp. 157–166, 2014.
  * Chiu et al. [2016] W.-Y. Chiu, G. G. Yen, and T.-K. Juan, “Minimum manhattan distance approach to multiple criteria decision making in multiobjective optimization problems,” _IEEE Transactions on Evolutionary Computation_ , vol. 20, no. 6, pp. 972–985, 2016.
  * Zhang et al. [2023c] H. Zhang, Q. Chen, B. Xue, W. Banzhaf, and M. Zhang, “A semantic-based hoist mutation operator for evolutionary feature construction in regression,” _IEEE Transactions on Evolutionary Computation_ , 2023.
  * Cava et al. [2021] W. L. Cava, P. Orzechowski, B. Burlacu, F. O. de Franca, M. Virgolin, Y. Jin, M. Kommenda, and J. H. Moore, “Contemporary symbolic regression methods and their relative performance,” in _Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 1)_ , 2021\.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin, “XGBoost: A scalable tree boosting system,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2016, pp. 785–794.
  * Ke et al. [2017] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, “Lightgbm: A highly efficient gradient boosting decision tree,” _Advances in Neural Information Processing Systems_ , vol. 30, pp. 3146–3154, 2017.
  * Burlacu et al. [2020] B. Burlacu, G. Kronberger, and M. Kommenda, “Operon c++ an efficient genetic programming framework for symbolic regression,” in _Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion_ , 2020, pp. 1562–1570.
  * Wei et al. [2023] H. Wei, H. Zhuang, R. Xie, L. Feng, G. Niu, B. An, and Y. Li, “Mitigating memorization of noisy labels by clipping the model prediction,” in _International Conference on Machine Learning_. PMLR, 2023, pp. 36 868–36 886.
  * Smith-Miles and Geng [2020] K. Smith-Miles and X. Geng, “Revisiting facial age estimation with new insights from instance space analysis,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 5, pp. 2689–2697, 2020.
  * Wever et al. [2021] M. Wever, A. Tornede, F. Mohr, and E. Hüllermeier, “Automl for multi-label classification: Overview and empirical evaluation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 43, no. 9, pp. 3037–3054, 2021.
  * Zimmer et al. [2021] L. Zimmer, M. Lindauer, and F. Hutter, “Auto-pytorch: Multi-fidelity metalearning for efficient and robust autodl,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 43, no. 9, pp. 3079–3090, 2021\.
  * Mundhenk et al. [2021] T. Mundhenk, M. Landajuela, R. Glatt, C. P. Santiago, B. K. Petersen _et al._ , “Symbolic regression via deep reinforcement learning enhanced genetic programming seeding,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 24 912–24 923, 2021.
  * Kamienny et al. [2022] P.-A. Kamienny, S. d’Ascoli, G. Lample, and F. Charton, “End-to-end symbolic regression with transformers,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 10 269–10 281, 2022.
  * Wang and Mao [2021] Z. Wang and Y. Mao, “On the generalization of models trained with sgd: Information-theoretic bounds and implications,” in _International Conference on Learning Representations_ , 2021.
  * Zhang et al. [2021b] H. Zhang, A. Zhou, and H. Zhang, “An evolutionary forest for regression,” _IEEE Transactions on Evolutionary Computation_ , vol. 26, no. 4, pp. 735–749, 2021.

|  Hengzhe Zhang (M’23) is a PhD student in computer science at Victoria
University of Wellington, New Zealand. He received the B.Sc. degree in
software engineering from Xiangtan University, Hunan, China, in 2019, and the
M.Sc. degree in computer science from East China Normal University, Shanghai,
China, in 2022. His current research interests include symbolic regression,
genetic programming, evolution computation, and statistical machine learning.
---|---
|  Qi Chen (M’14) received the B.E. degree in Automation from the University
of South China, Hunan, China in 2005 and the M.E. degree in Software
Engineering from Beijing Institute of Technology, Beijing, China in 2007, and
the PhD degree in computer science in 2018 at Victoria University of
Wellington, New Zealand. Currently, she is a Senior Lecturer in Artificial
Intelligence in School of Engineering and Computer Science at VUW. Qi’s
research interests including machine learning, evolutionary computation,
feature selection, feature construction, transfer learning, domain adaptation
and statistical learning theory. She serves as a reviewer of international
conferences, including AAAI and IJCAI, and international journals, including
IEEE Transactions on Evolutionary Computation and IEEE Transactions on
Cybernetics.
---|---
|  Bing Xue (M’10-SM’21) received the Ph.D. degree in computer science from
the Victoria University of Wellington, Wellington, New Zealand, in 2014. She
is currently a Professor of Artificial Intelligence, Deputy Director of Centre
for Data Science and Artificial Intelligence, and Deputy Head of School in the
School of Engineering and Computer Science at Victoria University of
Wellington. She has over 300 papers published in fully refereed international
journals and conferences and her research focuses mainly on evolutionary
computation and machine learning. Dr Xue is currently the Chair of IEEE CIS
Evolutionary Computation Technical Committee and Editor of IEEE CIS
Newsletter. She has also served as associate editor of several international
journals, such as IEEE CIM, IEEE TEVC and ACM TELO. She is also a Fellow of
Engineering New Zealand.
---|---
|  Mengjie Zhang (M’04-SM’10-F’19) received the Ph.D. degree in computer
science from RMIT University, Melbourne, VIC, Australia, in 2000. He is
currently a Professor of Computer Science, the Director of Centre for Data
Science and Artificial Intelligence, Victoria University of Wellington, New
Zealand. His current research interests include genetic programming, image
analysis, feature selection and reduction, job-shop scheduling, and
evolutionary deep learning and transfer learning. He has published over 800
research papers in refereed international journals and conferences. He is a
Fellow of the Royal Society of New Zealand, a Fellow of Engineering New
Zealand, a Fellow of IEEE, and an IEEE Distinguished Lecturer.
---|---
|  Wolfgang Banzhaf (M’20) received the Dr.rer.nat (Ph.D.) degree from the
Department of Physics, Technische Hochschule Karlsruhe (currently, Karlsruhe
Institute of Technology), Karlsruhe, Germany, in 1985. He was the University
Research Professor with the Department of Computer Science, Memorial
University of Newfoundland, St. John’s, NL, Canada, where he served as the
Head of Department from 2003 to 2009 and from 2012 to 2016. He is the John R.
Koza Chair of Genetic Programming with the Department of Computer Science and
Engineering and a member of the BEACON Center for the Study of Evolution in
Action, Michigan State University, East Lansing, MI, USA. Studies of
selforganization and the field of Artificial Life are also of very much
interest to him. He has become more involved with network research as it
applies to natural and man-made systems. His research interests are in the
field of bioinspired computing, notably evolutionary computation, and complex
adaptive systems.
---|---

##  Appendix A Empirically Verify Sharpness Estimation

The theorem in Section 3.4 shows the equivalence between adding noise to input
layers and adding noise to weights. To verify this, we tested it on a neural
network with one hidden layer of 2000 neurons. The network is randomly
initialized, where weight $w$ is drawn from a normal distribution
$\mathcal{N}(0,1)$, and noise $\epsilon$ and $\epsilon^{\prime}$ are drawn
from normal distributions $\mathcal{N}(0,w_{k}^{2})$ and
$\mathcal{N}(0,x_{i}^{2})$, respectively, and added to the weights $w_{k}$ and
inputs of each layer $x_{i}$, i.e., $w_{k}^{\prime}=w_{k}+\epsilon$ or
$x_{i}^{\prime}=x_{i}+\epsilon^{\prime}$. To ensure reliable results,
experiments are executed 100 times, each round with 100 randomly generated
inputs with 2000 dimensions. The results are presented in Figure 20. They show
that the distribution of outputs for weight perturbation and input
perturbation is similar, thus empirically verifying the theorem in Section
3.4.

Figure 20: Distribution of outputs for different perturbation strategies.

##  Appendix B Hyperparameter Analysis

For SAM-GP, two hyperparameters need to be considered: the first one is the
number of perturbations to estimate the sharpness, and the second one is the
magnitude of perturbation added to the semantics of subtrees. This section
performs a hyperparameter analysis to study the impact of these
hyperparameters on the results.

###  B.1 Analysis of the Number of Perturbation Iterations

As illustrated in Section 5.3, the perturbation process incurs additional
computational costs. In this section, we investigate the impact of the number
of perturbations on the final performance. We explored four different
perturbation iterations, namely $\\{1,5,10,20\\}$. The statistical pairwise
comparison of test $R^{2}$ scores for different numbers of perturbation
iterations is presented in Table VI. The results indicate that using 10
iterations to estimate variance is sufficient to achieve good test $R^{2}$
scores, and its training time falls within acceptable ranges, as shown in
Figure 21. In the field of deep learning, it has been observed that increasing
the number of iterations to solve the sharpness estimation problem could lead
to changes in the optimal perturbation magnitude, but it did not result in
improved test performance [34]. In the neural network domain, Gaussian
perturbations with approximately 3 to 10 iterations are adequate to yield good
results [62]. This range changes to 10 to 20 in GP-based feature construction.

Figure 21: Distributions of training time for different numbers of perturbation iterations. TABLE VI: Statistical comparison of test $R^{2}$ scores using different numbers of perturbation iterations. | 5 | 10 | 20
---|---|---|---
1 | 4(+)/33($\sim$)/21(-) | 5(+)/27($\sim$)/26(-) | 7(+)/25($\sim$)/26(-)
5 | — | 2(+)/47($\sim$)/9(-) | 5(+)/43($\sim$)/10(-)
10 | — | — | 3(+)/54($\sim$)/1(-)

###  B.2 Analysis of the Magnitude of Perturbations

In the proposed method, the magnitude of perturbation is another
hyperparameter. In this section, we conduct comparative experiments with
standard deviations of $\\{0.1,0.3,0.5\\}$ to investigate the impact of the
magnitude of perturbation on controlling overfitting. Table VII presents
pairwise statistical comparisons of test $R^{2}$ scores regarding the
magnitude of perturbations. From the experimental results, it is evident that
selecting an appropriate magnitude for the perturbations is crucial for SAM.
This finding is consistent with the application of SAM techniques in deep
learning algorithms [34]. Increasing the magnitude from 0.3 to 0.5 leads to a
deterioration in results on 15 datasets and improves results on only 6
datasets. Similarly, decreasing the magnitude from 0.3 to 0.1 worsens outcomes
on 24 datasets and improves results on only 5 datasets. Overall, the
experimental results suggest that 0.3 is a suitable starting point if not
using cross-validation to select the optimal magnitude of perturbation.

TABLE VII: Statistical comparison of test $R^{2}$ scores using different magnitudes of perturbations. | 0.1 | 0.5
---|---|---
0.3 | 24(+)/29($\sim$)/5(-) | 15(+)/37($\sim$)/6(-)
0.1 | — | 13(+)/21($\sim$)/24(-)

##  Appendix C Comparison of Ensemble Size

For the ensemble version of SAM-GP, the ensemble size is a parameter worth
considering. This paper follows the convention of using 100 models in the
final ensemble [63]. In this section, we present the changes in test $R^{2}$
scores with varying ensemble sizes, as shown in Table VIII. The experimental
results indicate that using an ensemble size of 100 yields the best results,
but the difference between an ensemble size of 30 and 100 is not very
significant. In fact, with the greedy ensemble selection technique [9], it is
possible to bridge the gap between an ensemble size of 30 and 100, which
merits further exploration.

TABLE VIII: Statistical comparison of test $R^{2}$ scores using different ensemble sizes. | 10 | 30 | 100
---|---|---|---
5 | 0(+)/58($\sim$)/0(-) | 0(+)/58($\sim$)/0(-) | 0(+)/53($\sim$)/5(-)
10 | — | 0(+)/58($\sim$)/0(-) | 0(+)/56($\sim$)/2(-)
30 | — | — | 0(+)/58($\sim$)/0(-)

##  Appendix D Further Analysis of Sharpness Estimation

In this paper, several default settings are directly drawn from the deep
learning community, such as adaptive noise scaling [37] and the 1-sharpness
estimation strategy [34]. However, although they have a strong theoretical
foundation, as shown in Section 3.4, they are not the only choices for
sharpness estimation. This section provides a more in-depth analysis of these
strategies.

###  D.1 Comparison of Noise Types

In this paper, we use Gaussian noise $\mathcal{N}(0,1)$ to perturb GP trees.
Gaussian noise is chosen due to its simplicity for theoretical analysis and
its effectiveness in controlling overfitting in neural networks. However, in
practice, any type of noise can be used to perturb GP trees. In this section,
we consider three types of noise:

  * •

Uniform noise: Noise is sampled from a uniform distribution
$\mathcal{U}(-1,1)$.

  * •

Laplacian noise: Noise is sampled from a Laplacian distribution
$\mathcal{L}(0,1)$. The probability density distribution $f(\epsilon)$ is
given by:

$f(\epsilon)=\frac{1}{2}\exp\left(-|\epsilon|\right).$ (19)

  * •

Ensemble noise: Ensemble noise consists of three types of noise: Gaussian
noise, uniform noise, and Laplacian noise. In each round of perturbation, one
of these three noises is randomly selected using a random number
$\xi\in\mathcal{U}(0,1)$,

$\operatorname{Noise}=\begin{cases}\mathcal{N}(0,1)&\text{if
}0\leq\xi<\frac{1}{3}\\\ \mathcal{U}(-1,1)&\text{if
}\frac{1}{3}\leq\xi<\frac{2}{3}\\\ \mathcal{L}(0,1)&\text{if
}\frac{2}{3}\leq\xi\leq 1\\\ \end{cases}$ (20)

Then, noise values are sampled from the corresponding distribution.

The comparison of $R^{2}$ scores for these three types of noise is presented
in Table IX. The experimental results show that there is no significant
difference between using Gaussian noise and ensemble noise to perturb inputs.
When using Laplacian noise, there is a decrease in test $R^{2}$ scores, as
seen in the results. As for the uniform distribution, it has significant
advantages on 6 datasets and significantly worse performance on 19 datasets.
Therefore, in the real world, it is recommended to use cross-validation to
determine the noise type, but the most widely used Gaussian noise is a good
default type of noise to perturb semantics.

TABLE IX: Comparison of test $R^{2}$ scores using different selection operators. | Uniform | Laplace | Ensemble
---|---|---|---
Normal | 19(+)/33($\sim$)/6(-) | 10(+)/45($\sim$)/3(-) | 3(+)/54($\sim$)/1(-)
Uniform | — | 11(+)/24($\sim$)/23(-) | 8(+)/31($\sim$)/19(-)
Laplace | — | — | 2(+)/49($\sim$)/7(-)

###  D.2 Comaprison of M-Sharpness

In this paper, 1-sharpness is utilized for calculating sharpness due to its
widespread adoption in the deep learning community, and it offers better
regularization for inducing sparse models [34]. In practice, M can be assigned
any arbitrary value or even set equal to the number of samples. Therefore,
this section provides an empirical comparison and considers three settings:

  * •

$M=4$: In this method, 4 samples are randomly grouped as a mini-batch, and the
maximum sharpness is calculated separately for each mini-batch.

  * •

$M=N$: This corresponds to the traditional n-sharpness approach. Sharpness is
calculated on the entire batch, and the sharpness value in the iteration
corresponding to the maximum sharpness is selected as the final sharpness.
Formally, $S_{\text{n-SAM}}$ is defined as:

$S_{\text{n-SAM}}=\max_{k=1,\dots,K}\frac{1}{N}\sum_{i=1,\dots,N}(\tilde{Y_{i}^{k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}.$
(21)

  * •

Gaussian Model Perturbation (GMP) [62]: GMP uses Gaussian noise instead of the
gradient method to perturb neural networks to estimate sharpness. With $K$
rounds of perturbations $Y^{k}$, where $k\in[1,K]$, GMP optimizes the mean
perturbed loss instead of the maximum perturbed loss. Thus, the sharpness
$S_{\text{GMP}}$ is defined as:

$S_{\text{GMP}}=\frac{1}{N}\frac{1}{K}\sum_{i=1,\dots,N}\sum_{k=1,\dots,K}(\tilde{Y_{i}^{k}}-Y_{i})^{2}-(\hat{Y_{i}}-Y_{i})^{2}$
(22)

The experimental results are presented in Table X. The experimental results
demonstrate that using 1-sharpness yields the highest test $R^{2}$ scores,
followed by 4-SAM. Importantly, 1-SAM significantly outperforms n-SAM on 22
datasets while underperforming on 5 datasets. Thus, in conclusion, 1-SAM is
recommended in sharpness-aware minimization for evolutionary feature
construction compared to n-SAM, which is consistent with the findings in deep
learning [34]. The experimental results about GMP are interesting, as 1-SAM
significantly outperforms n-SAM on 23 datasets while underperforming on 4
datasets. These results demonstrate that in the situation of multi-round
random perturbation for estimating sharpness, taking the maximum sharpness
over multiple rounds is better than taking the average sharpness over these
rounds. In the future, this approach is worth validating in neural networks.

TABLE X: Comparison of test $R^{2}$ scores using different batch sizes for sharpness estimation. | 4-SAM | n-SAM | GMP
---|---|---|---
1-SAM | 6(+)/50($\sim$)/2(-) | 22(+)/31($\sim$)/5(-) | 23(+)/31($\sim$)/4(-)
4-SAM | — | 17(+)/38($\sim$)/3(-) | 19(+)/35($\sim$)/4(-)
n-SAM | — | — | 2(+)/56($\sim$)/0(-)

###  D.3 Effectiveness of Adaptive Sharpness

Adaptive sharpness optimization, based on the scale of weights, was initially
proposed in the context of deep learning and has demonstrated superior
performance compared to standard sharpness optimization [37]. In this paper,
we employ the adaptive sharpness strategy because it is supported by a solid
theoretical foundation, as shown in Section 3.4. However, it is still possible
to consider two alternative strategies:

  * •

Batch Adaptive: This strategy scales the noise based on the standard deviation
of a batch of inputs rather than scaling each input separately.

$\tilde{\psi}_{a}(X)=\psi_{a}(X)+\mathcal{N}(0,\sigma_{\psi_{a}}(X))$ (23)

  * •

Without Adaptive (W/O): This strategy does not utilize adaptive techniques. In
this case, random values are directly added to the inputs without any scaling.
In other words, noise is synthesized by the following equation:

$\tilde{\psi}_{a}(X)=\psi_{a}(X)+\mathcal{N}(0,\sigma^{2})$ (24)

The experimental results are presented in Table XI. The results indicate that
instance-wise adaptive sharpness estimation achieves significantly better
performance compared to adaptive sharpness based on the standard deviation of
the entire dataset or simply ignoring adaptiveness. Thus, it is recommended to
use the instance-wise adaptive sharpness strategy in sharpness-aware
minimization of evolutionary feature construction. It is worth noting that the
data are standardized before training. In cases where data are not
standardized, the adaptive sharpness strategy may be even more necessary to
achieve good results.

TABLE XI: Comparison of test $R^{2}$ scores using different adaptive sharpness strategies. | Batch | W/O
---|---|---
Instance | 12(+)/42($\sim$)/4(-) | 12(+)/41($\sim$)/5(-)
Batch | — | 1(+)/56($\sim$)/1(-)
TABLE XII: The Hyperparameter Spaces of ML Methods used in Grid Search [52]. Method |  Hyperparameters
---|---
AdaBoost |  {‘learning_rate’: (0.01, 0.1, 1.0, 10.0), ‘n_estimators’: (10, 100, 1000)}
KernelRidge |  {‘kernel’: (‘linear’, ‘poly’, ‘rbf’, ‘sigmoid’), ‘alpha’: (0.0001, 0.01, 0.1, 1), ‘gamma’: (0.01, 0.1, 1, 10)}
LGBM |  {‘n_estimators’: (10, 50, 100, 250, 500, 1000), ‘learning_rate’: (0.0001, 0.01, 0.05, 0.1, 0.2), ‘subsample’: (0.5, 0.75, 1), ‘boosting_type’: (‘gbdt’, ‘dart’, ‘goss’)}
MLP |  {‘activation’: (‘logistic’, ‘tanh’, ‘relu’), ‘solver’: (‘lbfgs’, ‘adam’, ‘sgd’), ‘learning_rate’: (‘constant’, ‘invscaling’, ‘adaptive’)}
RandomForest |  {‘n_estimators’: (10, 100, 1000), ‘min_weight_fraction_leaf’: (0.0, 0.25, 0.5), ‘max_features’: (‘sqrt’, ‘log2’, None)}
LR (SGD) |  {‘alpha’: (1e-06, 0.0001, 0.01, 1), ‘penalty’: (‘l2’, ‘l1’, ‘elasticnet’)}
XGB |  {‘n_estimators’: (10, 50, 100, 250, 500, 1000), ‘learning_rate’: (0.0001, 0.01, 0.05, 0.1, 0.2), ‘gamma’: (0, 0.1, 0.2, 0.3, 0.4), ‘subsample’: (0.5, 0.75, 1)}
Operon |  {‘population_size’: (500,), ‘pool_size’: (500,), ‘max_length’: (50,), ‘allowed_symbols’: (‘add,mul,aq,constant,variable’,), ‘local_iterations’: (5,), ‘offspring_generator’: (‘basic’,), ‘tournament_size’: (5,), ‘reinserter’: (‘keep-best’,), ‘max_evaluations’: (500000,)}
|  {‘population_size’: (500,), ‘pool_size’: (500,), ‘max_length’: (25,),
‘allowed_symbols’: (‘add,mul,aq,exp,log,sin,tanh,constant,variable’,),
‘local_iterations’: (5,), ‘offspring_generator’: (‘basic’,),
‘tournament_size’: (5,), ‘reinserter’: (‘keep-best’,), ‘max_evaluations’:
(500000,)}
|  {‘population_size’: (500,), ‘pool_size’: (500,), ‘max_length’: (25,),
‘allowed_symbols’: (‘add,mul,aq,constant,variable’,), ‘local_iterations’:
(5,), ‘offspring_generator’: (‘basic’,), ‘tournament_size’: (5,),
‘reinserter’: (‘keep-best’,), ‘max_evaluations’: (500000,)}
|  {‘population_size’: (100,), ‘pool_size’: (100,), ‘max_length’: (50,),
‘allowed_symbols’: (‘add,mul,aq,constant,variable’,), ‘local_iterations’:
(5,), ‘offspring_generator’: (‘basic’,), ‘tournament_size’: (3,),
‘reinserter’: (‘keep-best’,), ‘max_evaluations’: (500000,)}
|  {‘population_size’: (100,), ‘pool_size’: (100,), ‘max_length’: (25,),
‘allowed_symbols’: (‘add,mul,aq,exp,log,sin,tanh,constant,variable’,),
‘local_iterations’: (5,), ‘offspring_generator’: (‘basic’,),
‘tournament_size’: (3,), ‘reinserter’: (‘keep-best’,), ‘max_evaluations’:
(500000,)}
|  {‘population_size’: (100,), ‘pool_size’: (100,), ‘max_length’: (25,),
‘allowed_symbols’: (‘add,mul,aq,constant,variable’,), ‘local_iterations’:
(5,), ‘offspring_generator’: (‘basic’,), ‘tournament_size’: (3,),
‘reinserter’: (‘keep-best’,), ‘max_evaluations’: (500000,)}
FFX |  N/A
TABLE XIII: Detailed mean test $R^{2}$ of top-performing algorithms on all 58 PMLB datasets. Some extremely poor results with values lower than -1 have been replaced with ”NaN”. | SAM-EGP | SAM-GP | AdaBoost | RandomForest | XGB | MLP | LGBM
---|---|---|---|---|---|---|---
1027 | 0.854726 | 0.854561 | 0.787393 | 0.818687 | 0.766583 | 0.855499 | 0.797568
1028 | 0.308487 | 0.290274 | 0.291492 | 0.247797 | 0.072232 | 0.211904 | 0.173838
1029 | 0.532151 | 0.532176 | 0.416140 | 0.399828 | 0.292998 | 0.485666 | 0.351490
1030 | 0.342372 | 0.343076 | 0.218239 | 0.193312 | 0.029510 | 0.284981 | 0.195258
1089 | 0.710169 | 0.693212 | 0.663939 | 0.638205 | 0.494464 | 0.730050 | 0.490654
1096 | 0.782031 | 0.779123 | 0.654888 | 0.739644 | 0.631863 | 0.798300 | 0.590756
1191 | 0.168656 | 0.134398 | 0.180747 | 0.202281 | 0.042035 | -0.092548 | 0.065320
1193 | 0.523828 | 0.499217 | 0.507364 | 0.514330 | 0.396540 | 0.438033 | 0.451329
1196 | 0.406076 | 0.384779 | 0.344141 | 0.352092 | 0.191686 | 0.256396 | 0.264480
1199 | 0.390203 | 0.381523 | 0.318973 | 0.346348 | 0.227642 | 0.248171 | 0.210688
1201 | 0.004921 | -0.035013 | -0.118123 | -0.111686 | -0.360642 | -0.232617 | -0.005870
1203 | 0.503840 | 0.478036 | 0.442429 | 0.486452 | 0.355591 | 0.407642 | 0.387054
1595 | -0.031270 | -0.076824 | -0.108539 | -0.087267 | -0.316255 | -0.404611 | -0.004469
192 | 0.467214 | 0.401665 | 0.502347 | 0.505674 | 0.345272 | 0.560730 | 0.323541
195 | 0.782467 | 0.750412 | 0.823562 | 0.846289 | 0.789803 | 0.784240 | 0.763104
197 | 0.965190 | 0.956489 | 0.944252 | 0.938034 | 0.942572 | 0.527992 | 0.934283
201 | 0.750311 | 0.756006 | 0.760735 | 0.786726 | 0.749537 | 0.287320 | 0.736526
207 | 0.770959 | 0.749170 | 0.804841 | 0.842715 | 0.783954 | 0.794989 | 0.756793
210 | 0.820529 | 0.803585 | 0.735931 | 0.768603 | 0.768930 | 0.837155 | 0.669371
215 | 0.910450 | 0.889931 | 0.876456 | 0.893822 | 0.888945 | 0.840620 | 0.882717
218 | 0.431111 | 0.416437 | 0.407862 | 0.443194 | 0.352278 | 0.247238 | 0.271315
225 | 0.546320 | 0.504644 | 0.469742 | 0.492196 | 0.403452 | 0.186283 | 0.406376
227 | 0.952632 | 0.946648 | 0.932127 | 0.931276 | 0.931976 | 0.621185 | 0.921244
228 | 0.559855 | 0.488469 | 0.670801 | 0.671009 | 0.636026 | 0.678584 | 0.617867
229 | 0.841345 | 0.825342 | 0.795461 | 0.811418 | 0.800843 | 0.785734 | 0.783997
230 | 0.753547 | 0.790716 | 0.761560 | 0.772730 | 0.715084 | 0.794769 | 0.750850
294 | 0.766088 | 0.735205 | 0.742512 | 0.750332 | 0.662901 | 0.753038 | 0.704116
344 | 0.996610 | 0.995030 | 0.916482 | 0.937511 | 0.942911 | 0.924209 | 0.937431
4544 | 0.695669 | 0.673989 | 0.579617 | 0.583654 | 0.475654 | 0.380311 | 0.370879
485 | 0.466712 | 0.477391 | 0.461214 | 0.479891 | 0.355181 | 0.408427 | 0.375554
503 | 0.708665 | 0.705272 | 0.643500 | 0.664546 | 0.595600 | 0.695853 | 0.580727
505 | 0.996106 | 0.995856 | 0.979081 | 0.978621 | 0.982174 | 0.954379 | 0.970657
519 | 0.735646 | 0.722003 | 0.696015 | 0.651892 | 0.596941 | 0.732528 | 0.662674
522 | 0.161427 | 0.131922 | 0.190446 | 0.237460 | 0.095320 | 0.088440 | 0.087532
523 | 0.940811 | 0.938427 | 0.943045 | 0.945999 | 0.941825 | 0.936877 | 0.944799
527 | 0.952894 | 0.948503 | 0.717565 | 0.722000 | 0.672494 | NaN | 0.541663
529 | 0.774732 | 0.776661 | 0.495957 | 0.534237 | 0.478548 | 0.751388 | 0.502199
537 | 0.621150 | 0.607407 | 0.484567 | 0.507297 | 0.379761 | 0.581172 | 0.388682
542 | 0.167115 | 0.192564 | 0.431409 | 0.357981 | 0.145668 | 0.164405 | 0.034206
547 | 0.469542 | 0.463897 | 0.475554 | 0.486000 | 0.402307 | 0.418915 | 0.411154
556 | 0.866024 | 0.857930 | 0.844026 | 0.863950 | 0.821196 | 0.809733 | 0.787475
557 | 0.868733 | 0.867873 | 0.844251 | 0.891474 | 0.855027 | 0.848002 | 0.802575
560 | 0.979913 | 0.982663 | 0.942899 | 0.968217 | 0.967661 | 0.947857 | 0.952834
561 | 0.929144 | 0.924079 | 0.707248 | 0.776173 | 0.853230 | 0.851298 | 0.876368
562 | 0.952632 | 0.946648 | 0.932127 | 0.931276 | 0.931976 | 0.621185 | 0.921244
564 | 0.911372 | 0.900269 | 0.666915 | 0.665891 | 0.626455 | 0.731343 | 0.673700
573 | 0.965190 | 0.956489 | 0.944252 | 0.938034 | 0.942572 | 0.527992 | 0.934283
574 | 0.211307 | 0.176390 | 0.318173 | 0.324906 | 0.182810 | -0.007378 | 0.110001
659 | 0.622348 | 0.592591 | 0.511390 | 0.529914 | 0.254964 | 0.564836 | 0.307717
663 | 0.995971 | 0.993546 | 0.952503 | 0.976445 | 0.985812 | 0.991243 | 0.973437
665 | 0.305412 | 0.310292 | 0.163178 | 0.208788 | -0.098119 | 0.205274 | 0.072138
666 | 0.563532 | 0.542018 | 0.491946 | 0.508282 | 0.420062 | 0.539845 | 0.339060
678 | 0.275103 | 0.273002 | 0.044757 | 0.111072 | -0.176710 | 0.185851 | 0.020555
687 | 0.475916 | 0.400928 | 0.428099 | 0.468948 | 0.298187 | 0.433818 | 0.234628
690 | 0.959882 | 0.956200 | 0.951978 | 0.959725 | 0.958619 | 0.933151 | 0.959247
695 | 0.859767 | 0.855277 | 0.824942 | 0.835866 | 0.780648 | 0.785195 | 0.773405
706 | 0.577552 | 0.502538 | 0.503907 | 0.586631 | 0.474260 | 0.552942 | 0.389710
712 | 0.739643 | 0.737850 | 0.704985 | 0.691104 | 0.622566 | 0.750884 | 0.668710