# Principled Design of Translation, Scale, and Rotation Invariant Variation
Operators
for Metaheuristics

Ye Tian, Xingyi Zhang, Cheng He, Kay Chen Tan, and Yaochu Jin Y. Tian is with
the Key Laboratory of Intelligent Computing and Signal Processing of Ministry
of Education, Institutes of Physical Science and Information Technology, Anhui
University, Hefei 230601, China (email: field910921@gmail.com).X. Zhang is
with the Key Laboratory of Intelligent Computing and Signal Processing of
Ministry of Education, School of Computer Science and Technology, Anhui
University, Hefei 230601, China (email: xyzhanghust@gmail.com).C. He is with
the Guangdong Provincial Key Laboratory of Brain-inspired Intelligent
Computation, Department of Computer Science and Engineering, Southern
University of Science and Technology, Shenzhen 518055, China (email:
chenghehust@gmail.com).K. C. Tan is with the Department of Computing, The Hong
Kong Polytechnic University, Hong Kong SAR (email: kctan@polyu.edu.hk).Y. Jin
is with the Department of Computer Science, University of Surrey, Guildford,
Surrey, GU2 7XH, U.K. (email: yaochu.jin@surrey.ac.uk).

###### Abstract

In the past three decades, a large number of metaheuristics have been proposed
and shown high performance in solving complex optimization problems. While
most variation operators in existing metaheuristics are empirically designed,
this paper aims to design new operators automatically, which are expected to
be search space independent and thus exhibit robust performance on different
problems. For this purpose, this work first investigates the influence of
translation invariance, scale invariance, and rotation invariance on the
search behavior and performance of some representative operators. Then, we
deduce the generic form of translation, scale, and rotation invariant
operators. Afterwards, a principled approach is proposed for the automated
design of operators, which searches for high-performance operators based on
the deduced generic form. The experimental results demonstrate that the
operators generated by the proposed approach outperform state-of-the-art ones
on a variety of problems with complex landscapes and up to 1000 decision
variables.

###### Index Terms:

Variation operator, translation invariance, scale invariance, rotation
invariance, automated design.

##  I Introduction

Metaheuristics have shown effectiveness in solving the complex optimization
problems from various fields, such as manufacturing [1], scheduling [2],
bioinformatics [3], and economics [4]. In contrast to mathematical programming
methods that require an explicit objective function [5], metaheuristics
provide a high-level methodology solving problems in a black-box manner. So
far, there have been a variety of algorithms inspired by the biological
mechanisms in natural evolution, such as genetic algorithms (GA) [6],
differential evolution (DE) [7], evolution strategies [8], and evolutionary
programming [9]. Moreover, many swarm intelligence based algorithms have also
been proposed, including particle swarm optimization (PSO) [10], ant colony
optimization [11], artificial bee colony algorithms [12], among many others.

Existing metaheuristics exhibit quite different search behaviors and
optimization performance, which are mainly determined by the manually designed
variation operators [13, 14, 15], i.e., the strategies for receiving the
decision vectors of parents and generating the decision vectors of offspring
solutions. For example, the GA was designed according to the evolution theory
and law of inheritance, which generates offspring by the crossover between two
parents and the mutation on a single offspring solution. The crossover and
mutation operators provide a powerful exploration ability [16], making GA good
at handling multimodal landscapes [17]. The DE mutates each solution according
to the weighted difference between the other two solutions, which holds a good
performance on problems with complicated variable linkages [18]. The
covariance matrix adaptation evolution strategy (CMA-ES) generates new
solutions by sampling a multivariate normal distribution model adaptively
learned from the population, showing high performance on many real-world
applications [19]. Inspired by the choreography of bird flocking, the PSO
updates each particle according to its personal best particle and the global
best particle, which has a high speed of convergence [20].

Among the superiorities in existing variation operators, the independence of
search space is crucial to the robustness and generalization of
metaheuristics, since metaheuristics do not rely on specific characteristics
of problems. A search space independent operator holds the same performance on
a problem with arbitrary search space transformations, including translation
(i.e., $x^{\prime}=x+b$), scaling (i.e., $x^{\prime}=ax$), and rotation (i.e.,
$\mathbf{x}^{\prime}=\mathbf{x}M$). Existing studies have shed some light on
the sufficient conditions of achieving these invariance properties. For
instance, CMA-ES is scale invariant since its step-sizes are set
proportionally to the distance to the optimum found so far [21], and the
mutation operator of DE is rotation invariant since it is the weighted sum of
multiple parents [22]. Nevertheless, the mathematical definitions of all the
three properties have not been given, and the sufficient and necessary
condition of achieving them is still unknown. Therefore, it is difficult to
analyze whether an operator is translation, scale, and rotation invariant
theoretically, or to consider these properties in designing new operators
explicitly.

To address this issue, this work aims to deduce the generic form (i.e.,
sufficient and necessary condition) of translation, scale, and rotation
invariant operators, which can be used to judge the possession of invariance
properties and guide the design of new operators. Based on the deduced generic
form, this work proposes a principled approach for designing new operators
with invariance properties. In contrast to the parameter tuning of
metaheuristics [23, 24], the off-line recommendation of metaheuristics [25],
and the on-line combination of metaheuristics [26], the proposed approach is
not based on any existing metaheuristic but searches for totally new
operators. Moreover, the proposed approach does not utilize any existing
optimizer or classifier (e.g., F-Race for parameter tuning [27], artificial
neural network for metaheuristic recommendation [28], and sum-of-ranks
multiarmed bandit algorithm for metaheuristic combination [29]). By contrast,
it is a self-contained approach that can search for new operators by itself.
The main components of this work include the following three aspects:

  * •

Theoretical analysis: To illustrate the importance of translation invariance,
scale invariance, and rotation invariance, their effects on the search
behavior and performance of metaheuristics are investigated. Then, the
sufficient and necessary condition of achieving these properties is
mathematically derived, which reveals the generic form of search space
independent operators.

  * •

New approach: A principled approach to automated design of variation operators
is proposed, termed AutoV. Based on the deduced generic form of variation
operators, AutoV converts the search of high-performance operators into an
optimization problem, in which the decision variables are the parameters in
the operators. This way, AutoV can solve optimization problems without relying
on any existing operators.

  * •

Experimental study: The variation operator found by AutoV is embedded in a
simple evolutionary framework and compared to eight classical or state-of-the-
art metaheuristics. The experimental results show that the operator found by
AutoV can obtain the best results on various challenging benchmark problems;
in particular, it outperforms the winner of the CEC competition that contains
multiple operators with complex adaptation strategies. The results indicate
that AutoV has the potential to replace the laborious process of manual design
of new metaheuristics.

The rest of this paper is organized as follows. Section II analyzes the
effects of the three invariance properties, and Section III deduces their
sufficient and necessary condition. Section IV presents the proposed
principled approach, and Section V gives the experimental studies. Finally,
Section VI concludes this paper.

##  II Effects of Invariance Properties

This work focuses on the variation operators for the following continuous
optimization problem:

$\begin{aligned} \min&\ \ \ f(\mathbf{x})\\\ {\rm s.\,t.}&\ \ \
\mathbf{l}\leq\mathbf{x}\leq\mathbf{u}\end{aligned}\ ,$ (1)

where $\mathbf{x}=(x_{1},x_{2},\dots,x_{D})$ is a decision vector denoting a
solution for the problem, $\mathbf{l}=(l_{1},l_{2},\dots,l_{D})$ denotes the
lower bound, $\mathbf{u}=(u_{1},u_{2},\dots,u_{D})$ denotes the upper bound,
and $D$ is the number of decision variables. To solve such problems in a
black-box manner, a variety of variation operators have been proposed to
generate solutions without using any specific information except for the lower
and upper bounds. This section first introduces some representative variation
operators, then presents the definitions of translation invariance, scale
invariance, and rotation invariance and analyzes their effects by examples and
experimental studies.

###  II-A Variation Operators in Metaheuristics

An operator generally receives the decision vectors of one or more parents,
then outputs the decision vectors of one or more offspring solutions. For
example, the simulated binary crossover (SBX) operator [30] used in GA uses
two parents $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ to generate two offspring
solutions $\mathbf{o}_{1}$ and $\mathbf{o}_{2}$ each time:

$\left\\{\begin{aligned}
o_{1d}=0.5\left[(1+\beta)x_{1d}+(1-\beta)x_{2d}\right]\\\
o_{2d}=0.5\left[(1-\beta)x_{1d}+(1+\beta)x_{2d}\right]\\\ \end{aligned}\ ,\
1\leq d\leq D\right.,$ (2)

where $o_{1d}$ denotes the $d$-th variable of solution $\mathbf{o}_{1}$ and
$\beta$ is a random number obeying a special distribution. The mutation
operator of DE/rand/1/bin [7] uses three parents $\mathbf{x}_{1}$,
$\mathbf{x}_{2}$, and $\mathbf{x}_{3}$ to generate an offspring solution
$\mathbf{o}$ each time:

$o_{d}=x_{1d}+F\cdot(x_{2d}-x_{3d})\ ,\ 1\leq d\leq D\ ,$ (3)

where $F$ is a parameter controlling the amplification of the difference
between $\mathbf{x}_{2}$ and $\mathbf{x}_{3}$. In contrast to the random
parameter $\beta$ in SBX that varies on each dimension, the parameter $F$ in
DE is a predefined constant. The operator of CMA-ES [8] generates offspring by
sampling a multi-variate normal distribution:

$\mathbf{o}=\mathbf{x}_{m}+\sigma\cdot\mathcal{N}(\mathbf{0},C)\ ,$ (4)

where $\mathbf{x}_{m}$ is the weighted sum of all solutions, $\sigma$ is a
vector of iteratively updated step-sizes, and $C$ is a covariance matrix
updated according to the current population. In general, the initial $\sigma$
can be set to $0.6(\mathbf{u}-\mathbf{l})$ [31]. Similar to CMA-ES, the
operator of fast evolutionary programming (FEP) [9] generates offspring by
sampling a single-variate normal distribution:

$o_{d}=x_{d}+\eta_{d}\cdot\mathcal{N}(0,1)\ ,\ 1\leq d\leq D\ ,$ (5)

where $\mathbf{\eta}$ is a vector of self-adaptive standard deviations related
to each solution $\mathbf{x}$, whose elements can be initialized to 3 [9].

It can be found that these operators generate offspring by using distinct
formulas. Generally, an operator can be regarded as a function
$h(x_{1d},x_{2d},\dots)$ performed on each dimension $d$, where
$x_{1d},x_{2d},\dots$ contains the decision variables of parents, the lower
bound, and the upper bound. Note that other parameters (e.g., $\beta$ in SBX
and $C$ in CMA-ES) are ignored for simplicity. In the following, we
investigate how the invariance properties influence the search behavior and
performance of operators.

###  II-B Effects of Translation Invariance

According to [32], a variation operator $h(x_{1d},x_{2d},\dots)$ is invariant
to search space transformation $\mathcal{T}$ means that
$h(\mathcal{T}(x_{1d}),\mathcal{T}(x_{2d}),\dots)=\mathcal{T}(h(x_{1d},x_{2d},\dots))$.
A translation of the search space can be regarded as an addition of each
decision variable, i.e., $\mathcal{T}(x)=x+b$, hence the translation
invariance property can be defined as follows:

######  Definition 1 (Translation invariance).

A variation operator $h(x_{1d},x_{2d},\dots)$ is translation invariant if and
only if

$h(x_{1d}+b,x_{2d}+b,\dots)=h(x_{1d},x_{2d},\dots)+b$ (6)

holds for any real constant $b$.

It is not difficult to find that all the operators of SBX, DE, CMA-ES, and FEP
are translation invariant. In particular, the SBX operator described in (2)
can be rewritten as

$h_{sbx}(x_{1d},x_{2d})=x_{1d}+0.5(1\pm\beta)(x_{2d}-x_{1d})\ ,$ (7)

hence

$\displaystyle h_{sbx}(x_{1d}+b,x_{2d}+b)$ (8)
$\displaystyle\qquad=x_{1d}+b+0.5(1\pm\beta)(x_{2d}+b-x_{1d}-b)$
$\displaystyle\qquad=h_{sbx}(x_{1d},x_{2d})+b$

and the operator is translation invariant. By contrast, if the SBX operator is
modified to

$h_{sbx^{\prime}}(x_{1d},x_{2d})=0.1x_{1d}+0.5(1\pm\beta)(x_{2d}-x_{1d})\ ,$
(9)

then

$\displaystyle h_{sbx^{\prime}}(x_{1d}+b,x_{2d}+b)$ (10)
$\displaystyle\qquad=0.1x_{1d}+0.1b+0.5(1\pm\beta)(x_{2d}+b-x_{1d}-b)$
$\displaystyle\qquad\neq h_{sbx^{\prime}}(x_{1d},x_{2d})+b$

and the operator becomes not translation invariant. Obviously, the modified
SBX′ operator is likely to evolve the population toward the origin.

(a) SBX based GA (translation invariant) (b) SBX′ based GA

Figure 1: Convergence profiles of SBX and SBX′ based GAs with fixed random
number seed on $f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and without
translation.

To better illustrate this fact, Fig. 1 depicts the convergence profiles of SBX
and SBX′ based GAs on problem $f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and
without translation, where the random number seed is fixed and both the
population size and the number of generations are set to 10. It can be found
that the SBX based GA holds the same search behavior and can always converge
to the global optimums of the three problems, whereas the SBX′ based GA always
converges to the origin. In short, the SBX operator is translation invariant
but SBX′ is not. Furthermore, Table I lists the mean and standard deviation of
the minimum objective values obtained by SBX and SBX′ based GAs on six
benchmark problems [17], averaged over 30 runs. The six problems have the same
global optimum $(0,0,\dots,0)$, while the global optimum is changed to
$(-6,-6,\dots,-6)$ if the problems are translated by $x^{\prime}=x+6$. For the
original problems, the SBX′ based GA significantly outperforms the SBX based
GA since the former can quickly converge to the origin. While for the
translated problems, the performance of SBX′ based GA deteriorates
considerably since it cannot converge to the translated global optimum; by
contrast, the performance of SBX based GA keeps unchanged since the SBX
operator is translation invariant.

TABLE I: Minimum Objective Values Obtained by SBX and SBX′ Based GAs on Six Problems with and without Translation. Original problem | SBX based GA | SBX′ based GA
---|---|---
($x_{1},\dots,x_{30}\in[-10,10]$) | (translation invariant)
Schwefel’s Function 2.22 | 5.4310e+0 (2.13e+0) | 1.0965e-49 (1.29e-49)
Schwefel’s Function 2.21 | 4.5792e+0 (1.21e+0) | 8.3896e-49 (1.23e-48)
Quaric Function | 1.6425e+3 (1.03e+3) | 2.9477e-4 (2.66e-4)
Griewank Function | 4.7700e-1 (1.87e-1) | 0.0000e+0 (0.00e+0)
Ackley’s Function | 3.3728e+0 (3.75e-1) | 8.8818e-16 (0.00e+0)
Rastrigin’s Function | 4.2269e+1 (1.08e+1) | 0.0000e+0 (0.00e+0)
Translated problem | SBX based GA | SBX′ based GA
($x^{\prime}=x+6$) | (translation invariant)
Schwefel’s Function 2.22 | 5.2237e+0 (1.88e+0) | 1.0083e+12 (2.18e+12)
Schwefel’s Function 2.21 | 4.9401e+0 (1.24e+0) | 5.0000e+0 (0.00e+0)
Quaric Function | 1.9875e+3 (1.62e+3) | 2.3718e+5 (1.53e+4)
Griewank Function | 4.6346e-1 (1.17e-1) | 1.1523e+0 (4.76e-3)
Ackley’s Function | 3.5154e+0 (7.18e-1) | 1.2639e+1 (5.98e-4)
Rastrigin’s Function | 4.5695e+1 (1.01e+1) | 7.4895e+2 (1.83e-1)

###  II-C Effects of Scale Invariance

Similarly, a scaling of the search space can be regarded as a multiplication
of each decision variable, i.e., $\mathcal{T}(x)=ax$, hence the scale
invariance property can be defined as follows:

######  Definition 2 (Scale invariance).

A variation operator $h(x_{1d},x_{2d},\dots)$ is scale invariant if and only
if

$h(ax_{1d},ax_{2d},\dots)=a\cdot h(x_{1d},x_{2d},\dots)$ (11)

holds for any real constant $a$.

It can be found that the operators of SBX, DE, and CMA-ES are scale invariant
but the operator of FEP is not. Considering that the step-size $\sigma$ is
proportionally related to the decision space, the operator of CMA-ES described
in (4) can be rewritten as

$h_{cmaes}(x_{md},l_{d},u_{d})=x_{md}+(u_{d}-l_{d})\cdot\mathcal{N}(0,{\sigma^{\prime}_{d}}^{2}c)\
,$ (12)

where $\sigma^{\prime}_{d}$ is a parameter within $(0,1]$ and $c$ is related
to the covariance matrix $C$. Therefore,

$\displaystyle h_{cmaes}(ax_{md},al_{d},au_{d})$ (13)
$\displaystyle\qquad=ax_{md}+(au_{d}-al_{d})\cdot\mathcal{N}(0,{\sigma^{\prime}_{d}}^{2}c)$
$\displaystyle\qquad=a\cdot h_{cmaes}(x_{md},l_{d},u_{d})\ ,$

which means that the operator is scale invariant. On the other hand, the
operator of FEP described in (5) can be rewritten as

$h_{fep}(x_{d})=x_{d}+\mathcal{N}(0,\eta_{d}^{2})\ ,$ (14)

hence

$h_{fep}(ax_{d})=ax_{d}+\mathcal{N}(0,\eta_{d}^{2})\neq a\cdot h_{fep}(x_{d})\
,$ (15)

which means that the operator is not scale invariant.

(a) CMA-ES (scale invariant) (b) FEP

Figure 2: Convergence profiles of CMA-ES and FEP with fixed random number seed
on $f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and without scaling.

Fig. 2 plots the convergence profiles of CMA-ES and FEP on
$f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and without scaling. It can be seen
that the search behaviors of CMA-ES are the same on the three problems with
different scales, whereas the search behaviors of FEP are quite different.
Moreover, Table II presents the mean and standard deviation of the minimum
objective values obtained by CMA-ES and FEP on six benchmark problems.
Obviously, CMA-ES is competitive to FEP on the original problems, while CMA-ES
dramatically outperforms FEP on the scaled problems. This is because the
operator of CMA-ES is scale invariant and thus exhibits the same performance
on a problem with different scales; by contrast, the operator of FEP is not
scale invariant and thus is sensitive to the scales of problems.

TABLE II: Minimum Objective Values Obtained by CMA-ES and FEP on Six Problems with and without Scaling. Original problem | CMA-ES | FEP
---|---|---
($x_{1},\dots,x_{30}\in[-10,10]$) | (scale invariant)
Schwefel’s Function 2.22 | 1.8226e+1 (3.38e+0) | 1.6919e+1 (5.30e+0)
Schwefel’s Function 2.21 | 1.0000e+1 (0.00e+0) | 4.6087e+0 (2.29e-1)
Quaric Function | 2.5024e+1 (9.92e+0) | 1.4005e+3 (1.69e+3)
Griewank Function | 3.2242e-1 (5.74e-2) | 8.5491e-1 (1.18e-1)
Ackley’s Function | 2.7634e+0 (1.71e-1) | 5.7111e+0 (5.97e-1)
Rastrigin’s Function | 2.1628e+2 (4.87e+0) | 1.6904e+2 (2.52e+1)
Scaled problem | CMA-ES | FEP
($x^{\prime}=10x$) | (scale invariant)
Schwefel’s Function 2.22 | 1.9390e+1 (2.70e+0) | 4.6714e+10 (5.47e+10)
Schwefel’s Function 2.21 | 1.3107e+0 (1.36e-1) | 8.4873e+0 (5.56e-1)
Quaric Function | 2.6711e+1 (2.22e+1) | 3.8536e+5 (2.35e+4)
Griewank Function | 2.8246e-1 (9.96e-2) | 1.1561e+0 (8.63e-3)
Ackley’s Function | 3.0355e+0 (3.42e-1) | 1.3471e+1 (5.21e-1)
Rastrigin’s Function | 2.1651e+2 (1.40e+1) | 9.2797e+2 (3.64e+1)

###  II-D Effects of Rotation Invariance

A rotation of the search space can be regarded as a matrix multiplication of
each decision vector, i.e., $\mathcal{T}(\mathbf{x})=\mathbf{x}M$, hence the
rotation invariance property can be defined as follows:

######  Definition 3 (Rotation invariance).

A variation operator $h(\mathbf{x}_{1},\mathbf{x}_{2},\dots)$ is rotation
invariant if and only if

$h(\mathbf{x}_{1}M,\mathbf{x}_{2}M,\dots)=h(\mathbf{x}_{1},\mathbf{x}_{2},\dots)M$
(16)

holds for any orthogonal matrix $M$.

Note that here the decision variables on all dimensions
$\mathbf{x}_{1},\mathbf{x}_{2},\dots$ rather than those on a single dimension
$x_{1d},x_{2d},\dots$ should be considered. It can be deduced that the
mutation operator of DE is rotation invariant while the operators of SBX, CMA-
ES, and FEP are not. In particular, the mutation operator of DE described in
(3) can be rewritten as

$h_{de}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})=\mathbf{x}_{1}+F\cdot(\mathbf{x}_{2}-\mathbf{x}_{3})\
,$ (17)

hence

$\displaystyle h_{de}(\mathbf{x}_{1}M,\mathbf{x}_{2}M,\mathbf{x}_{3}M)=$
$\displaystyle\ \mathbf{x}_{1}M+F\cdot(\mathbf{x}_{2}M-\mathbf{x}_{3}M)$ (18)
$\displaystyle=$ $\displaystyle\
h_{de}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})M\ ,$

and the operator is rotation invariant. By contrast, since the SBX operator
described in (7) can be rewritten as

$h_{sbx}(\mathbf{x}_{1},\mathbf{x}_{2})=\mathbf{x}_{1}+(\mathbf{x}_{2}-\mathbf{x}_{1})B\
,$ (19)

where

$\footnotesize B=\left(\begin{array}[]{cccc}0.5(1\pm\beta_{1})&0&\cdots&0\\\
0&0.5(1\pm\beta_{2})&\cdots&0\\\ \cdots&\cdots&\cdots&\cdots\\\
0&0&0&0.5(1\pm\beta_{D})\\\ \end{array}\right)\normalsize\ .$ (20)

Therefore,

$h_{sbx}(\mathbf{x}_{1}M,\mathbf{x}_{2}M)=\mathbf{x}_{1}M+(\mathbf{x}_{2}-\mathbf{x}_{1})MB$
(21)

and

$h_{sbx}(\mathbf{x}_{1},\mathbf{x}_{2})M=\mathbf{x}_{1}M+(\mathbf{x}_{2}-\mathbf{x}_{1})BM\
.$ (22)

That is,
$h_{sbx}(\mathbf{x}_{1}M,\mathbf{x}_{2}M)=h_{sbx}(\mathbf{x}_{1},\mathbf{x}_{2})M$
holds only if $MB=BM$, i.e., $\beta_{1}=\beta_{2}=\dots=\beta_{D}$, which is
almost impossible since they are independent random numbers. Hence, the SBX
operator is not rotation invariant.

(a) DE (rotation invariant) (b) SBX based GA

Figure 3: Convergence profiles of DE and SBX based GA with fixed random number
seed on $f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and without rotation, where
$M$ is a randomly generated orthogonal matrix.

Fig. 3 shows the convergence profiles of the mutation based DE and SBX based
GA on $f(x_{1},x_{2})=x_{1}^{4}+x_{2}^{4}$ with and without rotation. As can
be seen, the search behavior of DE keeps unchanged on the two problems, where
the convergence profile is rotated together with the search space. On the
contrary, the search behavior of GA is unstable and the convergence profiles
are different on the two problems. As can be further observed from Table III,
the performance of DE is similar to and much better than GA on six benchmark
problems with and without rotation, respectively, which is consistent with the
facts that the mutation operator of DE is rotation invariant while the SBX
operator is not. In fact, the superiority of DE on problems with complicated
variable linkages [18] is mainly due to its rotation invariance property.

TABLE III: Minimum Objective Values Obtained by DE and SBX Based GA on Six Problems with and without Rotation. Original problem | DE | SBX based GA
---|---|---
($x_{1},\dots,x_{30}\in[-10,10]$) | (rotation invariant)
Schwefel’s Function 2.22 | 1.9799e+1 (2.08e+0) | 5.4310e+0 (2.13e+0)
Schwefel’s Function 2.21 | 1.5589e+0 (1.87e-1) | 4.5792e+0 (1.21e+0)
Quaric Function | 1.8303e+2 (9.67e+1) | 1.6425e+3 (1.03e+3)
Griewank Function | 4.9032e-1 (1.01e-1) | 4.7700e-1 (1.87e-1)
Ackley’s Function | 3.7861e+0 (2.96e-1) | 3.3728e+0 (3.75e-1)
Rastrigin’s Function | 2.4783e+2 (2.34e+1) | 4.2269e+1 (1.08e+1)
Rotated problem | DE | SBX based GA
($\mathbf{x}^{\prime}=\mathbf{x}M$) | (rotation invariant)
Schwefel’s Function 2.22 | 1.7401e+1 (3.38e+0) | 2.8069e+1 (4.32e+0)
Schwefel’s Function 2.21 | 1.5535e+0 (2.21e-1) | 2.0725e+0 (4.82e-1)
Quaric Function | 1.2403e+2 (7.80e+1) | 9.8479e+2 (1.24e+3)
Griewank Function | 4.6714e-1 (7.22e-2) | 4.9478e-1 (1.20e-1)
Ackley’s Function | 3.8477e+0 (2.24e-1) | 4.0805e+0 (5.49e-1)
Rastrigin’s Function | 2.4588e+2 (1.23e+1) | 2.5095e+2 (1.12e+1)

##  III Sufficient and Necessary Condition of Achieving Invariance Properties

It can be concluded from the above analysis that the three invariance
properties are critical to the robustness of operators, and we can judge
whether an operator possesses these properties according to their definitions.
However, it is still tricky to consider them in the design of new operators
explicitly. Therefore, this section deduces the generic form of translation,
scale, and rotation invariant operators by three steps.

###  III-A Theoretical Derivation

Firstly, let $\mathbf{w}=(w_{1},w_{2},\dots)=(x_{1d},x_{2d},\dots)$, then a
translation invariant operator should satisfy

$h(\mathbf{w}+b)=h(\mathbf{w})+b\ .$ (23)

To find the generic form of $h(\mathbf{w})$, it is first expanded in a first
order Taylor series at an arbitrary point $\mathbf{w}_{0}$:

$h(\mathbf{w})=h(\mathbf{w}_{0})+(\mathbf{w}-\mathbf{w}_{0})\nabla
h(\mathbf{w}_{0}^{\prime})\ ,$ (24)

where $\nabla=(\frac{\partial}{\partial w_{1}},\frac{\partial}{\partial
w_{2}},\dots)^{T}$ and $\mathbf{w}_{0}^{\prime}$ is an unknown point. Based on
(23), we have

$\displaystyle h(\mathbf{w}_{0})+(\mathbf{w}+b-\mathbf{w}_{0})\nabla
h(\mathbf{w}_{0}^{\prime})$ (25)
$\displaystyle\qquad=h(\mathbf{w}_{0})+(\mathbf{w}-\mathbf{w}_{0})\nabla
h(\mathbf{w}_{0}^{\prime\prime})+b\ ,$

where $\mathbf{w}_{0}^{\prime}$ and $\mathbf{w}_{0}^{\prime\prime}$ are
unknown points determined by $\mathbf{w}_{0}$. Since (25) holds for any $b$,
we only consider the components including $b$ in (25):

$bI\nabla h(\mathbf{w}_{0}^{\prime})=b\ ,$ (26)

where $I$ is a vector of ones. That is,

$I\nabla h(\mathbf{w}_{0}^{\prime})=1$ (27)

holds for any $\mathbf{w}_{0}^{\prime}$, which measn that

$\frac{\partial h}{\partial w_{1}}+\frac{\partial h}{\partial w_{2}}+\cdots=1\
.$ (28)

Obviously, (28) is a quasilinear first-order nonhomogeneous partial
differential equation [33]. Let $h=h(\mathbf{w})$ be the solution of (28)
determined by a function $g(\mathbf{w},h)=0$, then the following homogeneous
partial differential equation can be obtained:

$\frac{\partial g}{\partial h}+\frac{\partial g}{\partial
w_{1}}+\frac{\partial g}{\partial w_{2}}+\cdots=0\ ,$ (29)

and thus the following first integrals can be obtained:

$\left\\{\begin{aligned} h-w_{1}=&\ c_{1}\\\ w_{1}-w_{2}=&\ c_{2}\\\
w_{2}-w_{3}=&\ c_{3}\\\ \cdots\ \ \cdots\end{aligned}\right.,$ (30)

Therefore, the solution of (29) is

$g=g(h-w_{1},w_{1}-w_{2},w_{2}-w_{3},\dots)\ .$ (31)

Since $g=0$, there must exists a function $\psi$ such that

$h-w_{1}=\psi(w_{1}-w_{2},w_{2}-w_{3},\dots)\ ,$ (32)

which is equivalent to

$h(\mathbf{w})=w_{1}+\psi(w_{2}-w_{1},w_{3}-w_{2},\dots)\ .$ (33)

Moreover, it is obvious that (33) satisfies (23), hence (33) is a sufficient
and necessary condition of (23), and the following theorem can be given:

######  Theorem 1 (Translation invariant operator).

A continuously differentiable variation operator $h(x_{1d},x_{2d},\dots)$ is
translation invariant if and only if it has the following form:

$h(x_{1d},x_{2d},\dots)=x_{1d}+\psi(x_{2d}-x_{1d},x_{3d}-x_{2d},\dots)\ ,$
(34)

where $\psi$ can be any continuously differentiable function.

Secondly, a scale invariant operator should satisfy

$h(a\mathbf{w})=a\cdot h(\mathbf{w})\ ,$ (35)

according to (33), a translation and scale invariant operator should satisfy

$\displaystyle aw_{1}+\psi(aw_{2}-aw_{1},aw_{3}-aw_{2},\dots)$ (36)
$\displaystyle\qquad=aw_{1}+a\cdot\psi(aw_{2}-aw_{1},aw_{3}-aw_{2},\dots)\ .$

Let $\mathbf{v}=(v_{1},v_{2},\dots)=(w_{2}-w_{1},w_{3}-w_{2},\dots)$, we have

$\psi(a\mathbf{v})=a\cdot\psi(\mathbf{v}),$ (37)

and the first order Taylor series expansion at an arbitrary point
$\mathbf{v}_{0}$ is

$\displaystyle\psi(\mathbf{v}_{0})+(a\mathbf{v}-\mathbf{v}_{0})\nabla\psi(\mathbf{v}_{0}^{\prime})$
(38)
$\displaystyle\qquad=a\psi(\mathbf{v}_{0})+a(\mathbf{v}-\mathbf{v}_{0})\nabla\psi(\mathbf{v}_{0}^{\prime\prime})\
,$

where $\mathbf{v}_{0}^{\prime}$ and $\mathbf{v}_{0}^{\prime\prime}$ are
unknown points determined by $\mathbf{v}_{0}$. Since (38) holds for any $a$,
we only consider the components including $a$ in (38), that is,

$\psi(\mathbf{v}_{0})=\mathbf{v}_{0}\nabla\psi(\mathbf{v}_{0}^{\prime})$ (39)

must hold for any $\mathbf{v}_{0}$ and $\mathbf{v}_{0}^{\prime}$, which means
that

$v_{1}\frac{\partial\psi}{\partial v_{1}}+v_{2}\frac{\partial\psi}{\partial
v_{2}}+\cdots=\psi\ .$ (40)

Let $g(\mathbf{v},\psi)=0$, the following homogeneous partial differential
equation can be obtained:

$\psi\frac{\partial g}{\partial\psi}+v_{1}\frac{\partial g}{\partial
v_{1}}+v_{2}\frac{\partial g}{\partial v_{2}}+\cdots=0\ ,$ (41)

and thus the following first integrals can be obtained:

$\left\\{\begin{aligned} \ln\psi-\ln v_{1}=&\ c_{1}\\\ \ln v_{1}-\ln v_{2}=&\
c_{2}\\\ \ln v_{2}-\ln v_{3}=&\ c_{3}\\\ \cdots\ \
\cdots\end{aligned}\right.,$ (42)

Therefore, the solution of (41) is

$g=g(\ln\psi-\ln v_{1},\ln v_{1}-\ln v_{2},\ln v_{2}-\ln v_{3},\dots)\ .$ (43)

Since $g=0$, there must exists a function $\varphi$ such that

$\ln\psi-\ln v_{1}=\varphi(\ln v_{1}-\ln v_{2},\ln v_{2}-\ln v_{3},\dots)\ ,$
(44)

hence the generic form of $\psi(\mathbf{v})$ can be determined:

$\ln\psi(\mathbf{v})=\ln
v_{1}+\varphi\left(\ln\frac{v_{1}}{v_{2}},\ln\frac{v_{2}}{v_{3}},\dots\right)\
.$ (45)

Let
$\varphi(u_{1},u_{2},\dots)=\ln\phi(e^{\frac{1}{u_{1}}},e^{\frac{1}{u_{2}}},\dots)$,
then

$\psi(\mathbf{v})=v_{1}\phi\left(\frac{v_{2}}{v_{1}},\frac{v_{3}}{v_{2}},\dots\right)\
.$ (46)

According to (33), we have

$h(\mathbf{w})=w_{1}+(w_{2}-w_{1})\phi\left(\frac{w_{3}-w_{2}}{w_{2}-w_{1}},\frac{w_{4}-w_{3}}{w_{3}-w_{2}},\dots\right).$
(47)

Moreover, it is obvious that (47) satisfies both (23) and (35), hence (47) is
a sufficient and necessary condition of (23) and (35), and the following
theorem can be given:

######  Theorem 2 (Translation and scale invariant operator).

A continuously differentiable variation operator $h(x_{1d},x_{2d},\dots)$ is
translation and scale invariant if and only if it has the following form:

$\footnotesize
h(x_{1d},x_{2d},\dots)=x_{1d}+(x_{2d}-x_{1d})\phi\left(\frac{x_{3d}-x_{2d}}{x_{2d}-x_{1d}},\frac{x_{4d}-x_{3d}}{x_{3d}-x_{2d}},\dots\right),$
(48)

where $\phi$ can be any continuously differentiable function.

Thirdly, a rotation invariant operator should satisfy

$\footnotesize
h\left(\sum_{i=1}^{D}m_{id}x_{1i},\sum_{i=1}^{D}m_{id}x_{2i},\dots\right)=\sum_{i=1}^{D}m_{id}\cdot
h(x_{1i},x_{2i},\dots)\ ,$ (49)

where $m_{id}\in M$ and $d=1,\dots,D$. Let
$\phi(\mathbf{w})=\varphi(w_{1},\prod_{i=1}^{2}w_{i},\prod_{i=1}^{3}w_{i},\dots)$,
then (47) is equivalent to

$h(\mathbf{w})=w_{1}+(w_{2}-w_{1})\varphi\left(\frac{w_{3}-w_{2}}{w_{2}-w_{1}},\frac{w_{4}-w_{3}}{w_{2}-w_{1}},\dots\right)\
.$ (50)

According to (49), a translation, scale, and rotation invariant operator
should satisfy

$\displaystyle\sum_{i=1}^{D}m_{id}x_{1i}+\left(\sum_{i=1}^{D}m_{id}x_{2i}-\sum_{i=1}^{D}m_{id}x_{1i}\right)$
(51)
$\displaystyle\cdot\varphi\left(\frac{\sum_{i=1}^{D}m_{id}x_{3i}-\sum_{i=1}^{D}m_{id}x_{2i}}{\sum_{i=1}^{D}m_{id}x_{2i}-\sum_{i=1}^{D}m_{id}x_{1i}},\dots\right)$
$\displaystyle=\sum_{i=1}^{D}m_{id}\left[x_{1i}+(x_{2i}-x_{1i})\varphi\left(\frac{x_{3i}-x_{2i}}{x_{2i}-x_{1i}},\dots\right)\right]\
,$

which is equivalent to

$\displaystyle\sum_{i=1}^{D}m_{id}(x_{2i}-x_{1i})\cdot\varphi\left(\frac{\sum_{i=1}^{D}m_{id}(x_{3i}-x_{2i})}{\sum_{i=1}^{D}m_{id}(x_{2i}-x_{1i})},\dots\right)$
(52)
$\displaystyle=\sum_{i=1}^{D}m_{id}(x_{2i}-x_{1i})\varphi\left(\frac{m_{id}(x_{3i}-x_{2i})}{m_{id}(x_{2i}-x_{1i})},\dots\right)\
.$

Let
$\mathbf{u}_{i}=(u_{i1},u_{i2},\dots)=(m_{id}(x_{2i}-x_{1i}),m_{id}(x_{3i}-x_{2i}),\dots)$,
we have

$\sum_{i=1}^{D}u_{i1}\cdot\varphi\left(\frac{\sum_{i=1}^{D}u_{i2}}{\sum_{i=1}^{D}u_{i1}},\dots\right)=\sum_{i=1}^{D}u_{i1}\varphi\left(\frac{u_{i2}}{u_{i1}},\dots\right)\
,$ (53)

and the first order Taylor series expansion at an arbitrary point
$\mathbf{u}_{0}=(u_{01},u_{02},\dots)$ is

$\displaystyle\sum_{i=1}^{D}u_{i1}\cdot\left[\varphi(\mathbf{u}_{0})+\left(\frac{\sum_{i=1}^{D}u_{i2}}{\sum_{i=1}^{D}u_{i1}}-u_{01},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime})\right]$
(54)
$\displaystyle=u_{11}\cdot\left[\varphi(\mathbf{u}_{0})+\left(\frac{u_{12}}{u_{11}}-u_{01},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime\prime})\right]$
$\displaystyle+u_{21}\cdot\left[\varphi(\mathbf{u}_{0})+\left(\frac{u_{22}}{u_{21}}-u_{01},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime\prime\prime})\right]$
$\displaystyle+\cdots\ ,$

which can be simplified as

$\displaystyle\left(\sum_{i=1}^{D}u_{i2}-u_{01}\sum_{i=1}^{D}u_{i1},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime})\qquad$
(55)
$\displaystyle\qquad=\left(u_{12}-u_{01}u_{11},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime\prime})$
$\displaystyle\qquad+\left(u_{22}-u_{01}u_{21},\dots\right)\nabla\varphi(\mathbf{u}_{0}^{\prime\prime\prime})$
$\displaystyle\qquad+\cdots\ ,$

where $\mathbf{u}_{0}^{\prime},\mathbf{u}_{0}^{\prime\prime},\dots$ are
unknown points determined by $\mathbf{u}_{0}$. Since (55) holds for any
$u_{12},u_{22},\dots$, we have

$\nabla\varphi(\mathbf{u}_{0}^{\prime})=\nabla\varphi(\mathbf{u}_{0}^{\prime\prime})=\nabla\varphi(\mathbf{u}_{0}^{\prime\prime\prime})=\cdots=\mathbf{c}\
,$ (56)

and the form of $\varphi(\mathbf{u})$ can only be

$\varphi(\mathbf{u})=c_{0}+c_{1}u_{1}+c_{2}u_{2}+\cdots\ ,$ (57)

where $c_{0},c_{1},c_{2},\dots$ are constants. According to (50),

$h(\mathbf{w})=w_{1}+c_{0}(w_{2}-w_{1})+c_{1}(w_{3}-w_{2})+c_{2}(w_{4}-w_{3})+\cdots\
.$ (58)

Let

$\left\\{\begin{aligned} &\ 1-c_{0}=r_{1}\\\ &c_{0}-c_{1}=r_{2}\\\
&c_{1}-c_{2}=r_{3}\\\ &\cdots\ \ \cdots\\\ \end{aligned}\right.\ ,$ (59)

we have

$h(\mathbf{w})=r_{1}w_{1}+r_{2}w_{2}+r_{3}w_{3}+\cdots$ (60)

and $r_{1}+r_{2}+r_{3}+\cdots=1$. Moreover, it is obvious that (60) satisfies
(23), (35), and (49), hence (60) is a sufficient and necessary condition of
(23), (35), and (49), and the following theorem can be given:

######  Theorem 3 (Translation, scale, and rotation invariant operator).

A continuously differentiable variation operator $h(x_{1d},x_{2d},\dots)$ is
translation, scale, and rotation invariant if and only if it has the following
form:

$h(x_{1d},x_{2d},\dots)=r_{1}x_{1d}+r_{2}x_{2d}+r_{3}x_{3d}+\cdots\ ,$ (61)

where $r_{1},r_{2},r_{3},\dots$ can be any real constants satisfying
$r_{1}+r_{2}+r_{3}+\cdots=1$.

###  III-B Remarks

According to the above theorem, the following three corollaries can be given:

  1. 1.

An operator satisfying (61) is scale invariant.

  2. 2.

An operator satisfying (61) with $r_{1}+r_{2}+r_{3}+\cdots=1$ is scale and
translation invariant.

  3. 3.

An operator satisfying (61) with $r_{1},r_{2},r_{3},\dots$ being constants is
scale and rotation invariant.

The theoretical proofs of these corollaries are given in the following.

######  Corollary 1.

If a variation operator satisfying

$h(x_{1d},x_{2d},\dots)=r_{1}x_{1d}+r_{2}x_{2d}+r_{3}x_{3d}+\cdots\ ,$ (62)

then it is scale invariant.

###### Proof:

Since

$\displaystyle h(ax_{1d},ax_{2d},\dots)=$ $\displaystyle\
ar_{1}x_{1d}+ar_{2}x_{2d}+ar_{3}x_{3d}+\cdots$ (63) $\displaystyle=$
$\displaystyle\ ah(x_{1d},x_{2d},\dots)\ ,$

the operator is scale invariant. ∎

######  Corollary 2.

If a variation operator satisfying

$h(x_{1d},x_{2d},\dots)=r_{1}x_{1d}+r_{2}x_{2d}+r_{3}x_{3d}+\cdots$ (64)

with $r_{1}+r_{2}+r_{3}+\cdots=1$, then it is scale and translation invariant.

###### Proof:

Since

$\displaystyle h(x_{1d}+b,x_{2d}+b,\dots)$ (65)
$\displaystyle\quad=r_{1}(x_{1d}+b)+r_{2}(x_{2d}+b)+r_{3}(x_{3d}+b)+\cdots$
$\displaystyle\quad=r_{1}x_{1d}+r_{2}x_{2d}+r_{3}x_{3d}+\cdots+(r_{1}+r_{2}+r_{3}+\cdots)b$
$\displaystyle\quad=h(x_{1d},x_{2d},\dots)+b\ ,$

according to Corollary 1, the operator is scale and translation invariant. ∎

######  Corollary 3.

If a variation operator satisfying

$h(x_{1d},x_{2d},\dots)=r_{1}x_{1d}+r_{2}x_{2d}+r_{3}x_{3d}+\cdots$ (66)

with $r_{1},r_{2},r_{3},\dots$ being constants, then it is scale and rotation
invariant.

###### Proof:

Since $r_{1},r_{2},r_{3},\dots$ are constants, we have

$h(\mathbf{x}_{1},\mathbf{x}_{2},\dots)=r_{1}\mathbf{x}_{1}+r_{2}\mathbf{x}_{2}+r_{3}\mathbf{x}_{3}+\cdots\
.$ (67)

Therefore,

$\displaystyle h(\mathbf{x}_{1}M,\mathbf{x}_{2}M,\dots)$ (68)
$\displaystyle\qquad=r_{1}(\mathbf{x}_{1}M)+r_{2}(\mathbf{x}_{2}M)+r_{3}(\mathbf{x}_{3}M)+\cdots$
$\displaystyle\qquad=(r_{1}\mathbf{x}_{1})M+(r_{2}\mathbf{x}_{2})M+(r_{3}\mathbf{x}_{3})M+\cdots$
$\displaystyle\qquad=h(\mathbf{x}_{1},\mathbf{x}_{2},\dots)M\ ,$

according to Corollary 1, the operator is scale and rotation invariant. ∎

Generally, most existing operators including those in GA, DE, and CMA-ES meet
the second corollary and are scale and translation invariant. However, the
operators in GA and CMA-ES are not rotation invariant since the weights (i.e.,
$\beta$ in (2) and $\mathcal{N}(0,{\sigma^{\prime}_{d}}^{2}c)$ in (12)) vary
on different dimensions. By contrast, the mutation operator of DE is rotation
invariant since the weights (i.e., $F$ in (3)) keep unchanged on all
dimensions.

It should be noted that a rotation invariant operator does not necessarily
lead to a rotation invariant metaheuristic, and vice versa. For example, the
mutation operator of DE is rotation variant, but DE is not rotation invariant
due to the crossover operator with $CR<1$ [22]. The operator of CMA-ES is not
rotation invariant, but CMA-ES is rotation invariant due to the rotation angle
adaptive covariance matrix $C$ [8]. In fact, a rotation invariant operator may
not obtain good performance since offspring solutions can only be the linear
combinations of parents. In practice, the rotation invariance can be achieved
by rotating offspring solutions according to a covariance matrix learnt from
the population [8, 34].

##  IV Principled Approach for Designing Operators

The function of (61) reveals the generic form of translation, scale, and
rotation invariant operators, while the optimal values of the weights
$r_{1},r_{2},r_{3},\dots$ are not provided. Therefore, this section proposes a
principled approach for designing new operators, which searches for high-
performance operators by optimizing the weights.

###  IV-A Parameterization of Variation Operators

In order to introduce randomness, the proposed approach AutoV represents each
weight by an independent normal distribution, and it optimizes the mean and
variance of each normal distribution instead of the weight. Formally, the
proposed AutoV searches for the operators having the following form:

$\begin{aligned} h(x_{1d},\dots,x_{td})=&\ \sum_{i=1}^{t}r_{i}x_{id}\\\ {\rm
s.\,t.}\ \ \sum_{i=1}^{t}r_{i}=&\ 1,\
r_{i}\sim\mathcal{N}(\mu_{i},\sigma_{i}^{2})\end{aligned}\ ,$ (69)

where $\mu_{i}\in[-1,1]$ and $\sigma_{i}\in[0,1]$ are the parameters to be
optimized. Note that the value of $r_{i}$ keeps unchanged for $d=1,\dots,D$,
hence $r_{i}$ is still a constant for generating each offspring solution. In
the following, we demonstrate that an algorithm equipped with the above
operator can converge to the global optimum of continuous optimization
problems if given sufficient time.

######  Theorem 4.

Given a continuous function $f(\mathbf{x})$ and its global optimal solution
$\mathbf{x}^{\ast}$, for any $\epsilon>0$ and a metaheuristic based on an
elite strategy and a variation operator satisfying

$\begin{aligned} h(x_{1d},\dots,x_{td})=&\
(1-\sum_{i=2}^{t}r_{i})x_{1d}+\sum_{i=2}^{t}r_{i}x_{id}\\\ {\rm s.t.}\ \
r_{i+1}\sim&\ \mathcal{N}(\mu_{i},\sigma_{i}^{2}),\
i=1,\dots,t-1\end{aligned}\ ,$ (70)

we have

$\lim_{g\rightarrow\infty}P(|f(\mathbf{x}^{g})-f(\mathbf{x}^{\ast})|<\epsilon)=1\
,$ (71)

where $\mathbf{x}^{g}$ denotes the best solution in the population at the
$g$-th generation.

###### Proof:

Regarding the evolutionary procedure as a finite-dimension Markov chain [35],
let

$\begin{aligned} s_{1}:&\ |f(\mathbf{x}^{g})-f(\mathbf{x}^{\ast})|<\epsilon\\\
s_{2}:&\ |f(\mathbf{x}^{g})-f(\mathbf{x}^{\ast})|\geq\epsilon\end{aligned}\ ,$
(72)

then the one-step transition probabilities $P(s_{1}|s_{1})=1$ and
$P(s_{2}|s_{1})=0$, since the elite strategy always retains the best solution.
Now we focus on the transition probabilities $P(s_{1}|s_{2})$ and
$P(s_{2}|s_{2})$. Since $f$ is a continuous function, there must exist a
$\varepsilon>0$ such that $|f(\mathbf{x})-f(\mathbf{x}^{\ast})|<\epsilon$ when
$\|\mathbf{x}-\mathbf{x}^{\ast}\|_{\infty}<\varepsilon$, where
$\|\mathbf{x}-\mathbf{x}^{\ast}\|_{\infty}$ denotes the largest difference
between $\mathbf{x}$ and $\mathbf{x}^{\ast}$ over all the $D$ dimensions.
According to equation (70),

$\displaystyle h(x_{1d},\dots,x_{td})=$ $\displaystyle\
(1-\sum_{i=2}^{t}r_{i})x_{1d}+\sum_{i=2}^{t}r_{i}x_{id}$ (73) $\displaystyle=$
$\displaystyle\ x_{1d}+\sum_{i=2}^{t}r_{i}(x_{id}-x_{1d})$ $\displaystyle=$
$\displaystyle\ x_{1d}+\sum_{i=2}^{t}r_{i}^{\prime}\ ,$

where

$r^{\prime}_{i+1}\sim\mathcal{N}(\mu_{i},\sigma_{i}^{2}(x_{(i+1)d}-x_{1d})^{2})\
.$ (74)

Therefore,

$\displaystyle
P\left(\|h(\mathbf{x}_{1},\dots,\mathbf{x}_{t})-\mathbf{x}^{\ast}\|_{\infty}<\varepsilon\right)$
(75) $\displaystyle=$ $\displaystyle\
\prod_{d=1}^{D}P\left(|h(x_{1d},\dots,x_{td})-x_{d}^{\ast}|<\varepsilon\right)$
$\displaystyle=$ $\displaystyle\
\prod_{d=1}^{D}P\left(\left|x_{1d}+\sum_{i=2}^{t}r_{id}^{\prime}-x_{d}^{\ast}\right|<\varepsilon\right)$
$\displaystyle=$ $\displaystyle\
\prod_{d=1}^{D}P\left(x_{d}^{\ast}-x_{1d}-\varepsilon<\sum_{i=2}^{t}r_{id}^{\prime}<x_{d}^{\ast}-x_{1d}+\varepsilon\right)\
.$

Since all the $r_{id}^{\prime}$ are normally distributed,

$\displaystyle
P\left(\frac{x_{d}^{\ast}-x_{1d}-\varepsilon}{t-1}<r_{id}^{\prime}<\frac{x_{d}^{\ast}-x_{1d}+\varepsilon}{t-1}\right)$
(76) $\displaystyle=$ $\displaystyle\
\int_{\frac{x_{d}^{\ast}-x_{1d}-\varepsilon}{t-1}}^{\frac{x_{d}^{\ast}-x_{1d}+\varepsilon}{t-1}}\frac{e^{-\frac{c^{2}}{2\sigma_{i}^{2}(x_{(i+1)d}-x_{1d})^{2}}}}{\sigma_{i}(x_{(i+1)d}-x_{1d})\sqrt{2\pi}}dc>0\
,$

and thus

$P\left(\|h(\mathbf{x}_{1},\dots,\mathbf{x}_{t})-\mathbf{x}^{\ast}\|_{\infty}<\varepsilon\right)>0\
,$ (77)

which further indicates that $P(s_{1}|s_{2})>0$. Since
$P(s_{1}|s_{2})+P(s_{2}|s_{2})=1$, we have $P(s_{2}|s_{2})<1$. Therefore,

$\lim_{g\rightarrow\infty}\left[P(s_{2}|s_{2})\right]^{g}=0\ ,$ (78)

which means that the probability of
$s_{1}:|f(\mathbf{x}^{g})-f(\mathbf{x}^{\ast})|<\epsilon$ when
$g\rightarrow\infty$ is 1. ∎

It is worth noting that some existing operators include multiple parameter
sets to be selected with given probabilities (e.g., $\beta$ in SBX), which
provide complex search behaviors and can better balance between exploitation
and exploration. Hence, the ensemble of multiple parameter sets is also
adopted in AutoV. Moreover, the first weight $r_{1}$ can be omitted since it
can be directly set to $1-r_{2}-r_{3}-\cdots$. To summarize, an operator in
AutoV is determined by the following matrix

$\left(\begin{aligned}
\mu_{12},\sigma_{12},\mu_{13},\sigma_{13},\dots,\mu_{1t},\sigma_{1t},p_{1}\\\
\mu_{22},\sigma_{22},\mu_{23},\sigma_{23},\dots,\mu_{2t},\sigma_{2t},p_{2}\\\
\dots\ \ \dots\ \ \ \ \ \ \ \ \ \ \ \ \ \\\
\mu_{k2},\sigma_{k2},\mu_{k3},\sigma_{k3},\dots,\mu_{kt},\sigma_{kt},p_{k}\end{aligned}\right)\
,$ (79)

where $\mu_{ji}$ and $\sigma_{ji}$ denote the mean and variance of weight
$r_{i}$ in the $j$-th set, respectively, and $p_{j}$ denotes the probability
of selecting the $j$-th set. When generating a decision variable of an
offspring solution, the roulette-wheel selection is first used to select a
parameter set (i.e., one row of (79)) according to their probabilities. Then,
the decision variable is generated by sampling the given normal distribution.

In this way, the search of high-performance operators can be formulated as a
continuous optimization problem, whose decision vector contains all the
elements in (79) and the objective is the performance of the corresponding
operator. Generally, it is easy to measure the fitness by investigating the
performance of a metaheuristic equipped with the operator, but it is difficult
to optimize the decision vector since AutoV does not require any prior
knowledge, i.e., none of the existing optimizers are used to evolve a
population for solving the problem. To address this issue, the proposed AutoV
suggests a novel evolutionary procedure, in which the population can be
evolved by itself.

###  IV-B Procedure of the Principled Approach

Figure 4: Procedure of AutoV.

The procedure of the proposed AutoV is detailed in Fig. 4 and Algorithm 1.
Given a benchmark problem for performance measurement, AutoV first randomly
initializes a population $P$ and evaluates each solution in $P$, where each
solution denotes a parameter matrix defined in (79). At each generation, AutoV
selects a number of parents from $P$ via binary tournament selection, then
uses the parents to generate an offspring population by using the variation
operator parameterized by the best solution in $P$. Afterwards, the offspring
population is combined with $P$, and half the solutions with better fitness
survive for the next generation. The proposed AutoV repeats the above steps
until the termination criterion is fulfilled, and returns the solution with
the best fitness as the found variation operator.

Input: $f$ (a benchmark problem)

Output: $\mathbf{p}$ (the best solution)

1 $P\leftarrow$ Randomly initialize a population, where each solution is a
parameter matrix defined in (79);

2 for  _each $\mathbf{x}\in P$_ do

3 Evaluate the objective value of $\mathbf{x}$ by $Evaluation(\mathbf{x},f)$;

4

5while  _termination criterion is not fulfilled_ do

6 $P^{\prime}\leftarrow$ Select parents from $P$ via binary tournament
selection;

7 $\mathbf{p}\leftarrow$ The best solution in $P$;

8 $O\leftarrow$ Use the operator $\mathbf{p}$ to generate offspring based on
parents $P^{\prime}$;

9 for  _each $\mathbf{o}\in O$_ do

10 Evaluate the objective value of $\mathbf{o}$ by $Evaluation(\mathbf{o},f)$;

11 $P\leftarrow P\cup O$;

12 $P\leftarrow$ Half the solutions in $P$ with better fitness;

13

14$\mathbf{p}\leftarrow$ The best solution in $P$;

15 return $\mathbf{p}$;

Algorithm 1 Procedure of AutoV

Since the goal of AutoV is to create high-performance variation operators, it
uses the best operator found so far to evolve the population for finding a
better operator, and the better operator can then evolve the population for
finding a much better operator, where none of the existing metaheuristics are
used. For the fitness evaluation of each candidate operator, a simple
metaheuristic is established by adopting the candidate operator. As presented
in Algorithm 2, the best objective value found by the established
metaheuristic on the given benchmark problem is regarded as the fitness of the
candidate operator. Besides, we find that a candidate operator may be over-
optimized for the given benchmark problem, which means that it considers a
perfect solution for the given benchmark problem accidentally. Still, it
cannot evolve the population gradually. To solve this issue, AutoV executes
the established metaheuristic for multiple runs and uses the median value of
the found best objective values as the fitness.

Input: $\mathbf{p}$ (parameter matrix of an operator), $f$ (a benchmark
problem)

Output: $fit$ (fitness of $\mathbf{p}$)

1 $fit\leftarrow\emptyset$;

2 for  _$run=1$ to $maxRun$_ do

3 $P\leftarrow$ Randomly initialize a population, where each solution is a
decision vector for $f$;

4 for  _each $\mathbf{x}\in P$_ do

5 Evaluate the objective of $\mathbf{x}$ by $f(\mathbf{x})$;

6

7 while  _termination criterion is not fulfilled_ do

8 $P^{\prime}\leftarrow$ Select parents from $P$ via binary tournament
selection;

9 $O\leftarrow$ Use the operator $\mathbf{p}$ to generate offspring based on
parents $P^{\prime}$;

10 for  _each $\mathbf{o}\in O$_ do

11 Evaluate the objective of $\mathbf{o}$ by $f(\mathbf{o})$;

12 $P\leftarrow P\cup O$;

13 $P\leftarrow$ Half the solutions in $P$ with better fitness;

14

15 $\mathbf{x}\leftarrow$ The best solution in $P$;

16 $fit\leftarrow fit\cup\\{f(\mathbf{x})\\}$;

17

18$fit\leftarrow$ Median value of $fit$;

19 return $fit$;

Algorithm 2 $Evaluation(\mathbf{p},f)$

In the experiments, we consider the following five different sets of parents
as the input of (70) for generating one offspring solution:

$\left\\{\begin{aligned} h_{1}=&\ h(x_{1d},x_{2d})\\\ h_{2}=&\
h(x_{1d},l_{d},u_{d})\\\ h_{3}=&\ h(x_{1d},x_{2d},l_{d},u_{d})\\\ h_{4}=&\
h(x_{1d},x_{2d},x_{3d})\\\ h_{5}=&\ h(x_{1d},x_{2d},x_{3d},l_{d},u_{d})\\\
\end{aligned}\right.,$ (80)

where $x_{1d},x_{2d},x_{3d}$ denote the decision variables of three parents,
$l_{d}$ denotes the lower bound, and $u_{d}$ denotes the upper bound on the
$d$-th dimension.

##  V Experimental Studies

TABLE IV: Parameter Settings of the Compared Metaheuristics, Where $D$ is the Number of Decision Variables, $\mathbf{u}$ is the Upper Bound, $\mathbf{l}$ is the Lower Bound, and $n=100$ is the Population Size. Metaheuristic |  Parameter setting
---|---
GA [6] | crossover probability: 1,
(based on SBX [30] and |  mutation probability: $1/D$,
polynomial mutation [36]) | distribution index: 20
PSO [10] | inertia weight: 0.4
DE [7] | $CR=0.9$, $F=0.5$
(DE/rand/1/bin)
|  initial $\sigma=0.6(\mathbf{u}-\mathbf{l})$,
CMA-ES [8] |  $w_{i}=\log(\mu+0.5)-\log(i)$, $i\in[1,\mu]$,
$\mu=0.5n$, $c_{\sigma}=0.15$,
|  $d_{\sigma}=1.15$, $c_{c}=0.105$
FEP [9] |  initial $\eta=3$
CSO [37] | social factor: 0.1
(competitive swarm
optimizer)
SHADE [38] |  –
(parameter adaptative DE)
IMODE [39] | Minimum population size: 4
(winner of CEC’2020) | Ratio of archive size: 2.6

To verify the effectiveness of the proposed AutoV, the performance of the
operators found by AutoV is first studied. Then, the best operator is compared
with eight metaheuristics, where GA [6], PSO [10], DE [7], CMA-ES [8], and FEP
[9] are classical metaheuristics, CSO [37] is a competitive swarm optimizer
for large-scale optimization, and SHADE [38] and IMODE [39] are hybrid
metaheuristics based on multiple operators with parameter adaptation. Based on
the settings suggested in the original literature of the compared
metaheuristics, we finely tune their parameters for a relatively good
performance, where the detailed parameter settings are listed in Table IV.

###  V-A Comparison Between the Operators Found by AutoV

For each of the five functions given in (80), the proposed AutoV optimizes it
with a population size of 100 for 1000 generations. As for the fitness
evaluation of each candidate operator, the population size is set to 100, the
number of generations is set to 100, the number of parameter sets $k$ is set
to 10, and the number of runs $maxRun$ is set to 9. Besides, the Rastrigin’s
Function [17] is adopted as the benchmark problem for performance measurement.

To compare the performance of the found operators with different functions
$h_{1}$–$h_{5}$, they are embedded in the simple metaheuristic presented in
Algorithm 2 and tested on eight benchmark problems with 30 decision variables.
These benchmark problems have a variety of unimodal, multimodal, or flat
landscapes, whose definitions can be found in [17]. For all the
metaheuristics, the population size is set to 100 and the number of
generations is set to 100. Table V lists the minimum objective values found by
the five metaheuristics averaged over 30 runs, where the compared operators
exhibit similar performance and the operator $h_{3}$ has slightly better
overall performance than the others. It is worth noting that the operator
$h_{5}$ does not obtain the best overall performance, though the function
$h_{5}$ has more parents and is expected to perform better. This is because
the function $h_{5}$ contains more parameters to be optimized, which hinders
AutoV from finding high-performance operators.

TABLE V: Minimum Objective Values Obtained by the Operators with Different Functions $h_{1}$–$h_{5}$ Found by AutoV on Eight Benchmark Problems. Best Results are Highlighted. Problem | $h_{1}$ | $h_{2}$ | $h_{3}$ | $h_{4}$ | $h_{5}$
---|---|---|---|---|---
Schwefel’s | .07e-5 | 5.27e+0 | 2.07e-2 | 8.72e-3 | 1.02e-1
Function 2.22
Schwefel’s | 1.65e+1 | 1.19e+1 | 1.66e+0 | 8.75e+0 | 1.67e+0
Function 2.21
Quartic | 2.01e-2 | 1.32e-1 | 1.12e-2 | 9.02e-3 | 1.16e-2
Function
Generalized |  |  |  |  |
Griewank | 1.99e-3 | 2.93e+0 | 2.81e-2 | 4.46e-1 | 2.66e-1
Function |  |  |  |  |
Generalized |  |  |  |  |
Schwefel’s | -9.17e+3 | -1.04e+4 | -1.16e+4 | -5.70e+3 | -5.28e+3
Function 2.26 |  |  |  |  |
Ackley’s | 5.47e-4 | 4.77e+0 | 1.79e-2 | 8.94e-1 | 2.05e-1
Function
Rosenbrock’s | 1.20e+4 | 1.17e+5 | 9.09e+3 | 8.80e+3 | 3.82e+3
Function
Rastrigin’s | 3.76e+1 | 3.96e+1 | 3.88e+0 | 9.99e+0 | 7.86e+0
Function
(a)  Figure 5: Minimum objective values obtained by the operator with function $h_{3}$ and different numbers of parameter sets $k$ found by AutoV. TABLE VI: Minimum Objective Values Obtained by the Operator with Function $h_{3}$ and Different Benchmark Problems for Fitness Evaluation. Best Results are Highlighted. Problem | Optimized on | Optimized on | Optimized on | Optimized on
---|---|---|---|---
Schwefel’s | Schwefel’s | Ackley’s | Rastrigin’s
Function 2.22 | Function 2.21 | Function | Function
Schwefel’s | 2.46e-6 | 7.00e-2 | 4.91e-6 | 2.07e-2
Function 2.22
Schwefel’s | 1.46e+1 | 1.36e+0 | 1.21e+1 | 1.66e+0
Function 2.21
Quartic | 1.29e-2 | 1.23e-2 | 1.44e-2 | 1.12e-2
Function
Generalized | 6.89e-3 | 7.02e-2 | 1.48e-3 | 2.81e-2
Griewank
Function
Generalized | -9.35e+3 | -9.76e+3 | -9.55e+3 | -1.16e+4
Schwefel’s
Function 2.26
Ackley’s | 1.00e-4 | 5.77e-2 | 5.71e-5 | 1.79e-2
Function
Rosenbrock’s | 6.03e+3 | 3.89e+3 | 7.64e+3 | 9.09e+3
Function
Rastrigin’s | 2.36e+1 | 1.59e+1 | 2.91e+1 | 3.88e+0
Function
TABLE VII: Minimum Objective Values Obtained by Nine Metaheuristics on 13
Small-Scale Benchmark Problems with 30 Decision Variables. Best Results are
Highlighted.

Small-scale | GA [6] | PSO [10] | DE [7] | CMA-ES [8] | FEP [9] | CSO [37] | SHADE [38] | IMODE [39] | AutoV
---|---|---|---|---|---|---|---|---|---
problems [9]
$f_{1}$ | 1.4973e+1 $-$ | 3.5453e+3 $-$ | 5.0032e+3 $-$ | 1.2391e+2 $-$ | 9.2809e+3 $-$ | 1.2988e+3 $-$ | 2.2958e+1 $-$ | 2.9105e+0 $-$ | 8.3702e-1
(2.0053e+0) | (1.2242e+3) | (5.7205e+2) | (2.6004e+1) | (2.7503e+3) | (5.0772e+2) | (7.4054e+0) | (9.0907e-1) | (5.333e-1)
$f_{2}$ | 8.3310e-1 $-$ | 3.2170e+1 $-$ | 5.1595e+1 $-$ | 8.6103e+0 $-$ | 5.1291e+1 $-$ | 1.5504e+1 $-$ | 6.2312e+0 $-$ | 3.2305e-1 $-$ | 1.2081e-1
(1.7114e-1) | (6.6473e+0) | (1.5480e+1) | (3.3770e+0) | (1.0920e+1) | (3.1662e+0) | (1.1158e+0) | (1.0412e-1) | (2.556e-2)
$f_{3}$ | 1.0621e+4 $-$ | 7.8076e+3 $-$ | 6.4156e+3 $-$ | 5.0816e+3 $-$ | 2.1959e+4 $-$ | 1.7542e+3 $\approx$ | 1.6461e+3 $\approx$ | 6.6967e+2 $+$ | 2.0131e+3
(2.6732e+3) | (2.9985e+3) | (1.1614e+3) | (1.2254e+3) | (3.5707e+3) | (5.5302e+2) | (5.3479e+2) | (1.3529e+2) | (8.307e+2)
$f_{4}$ | 2.3222e+1 $-$ | 2.4005e+1 $-$ | 3.4503e+1 $-$ | 8.7765e+0 $-$ | 5.3464e+1 $-$ | 1.1521e+1 $-$ | 9.1224e+0 $-$ | 8.8346e+0 $-$ | 6.4179e+0
(3.6870e+0) | (3.0366e+0) | (4.1625e+0) | (1.5806e+0) | (9.3989e+0) | (2.0034e+0) | (1.4948e+0) | (2.0574e+0) | (1.953e+0)
$f_{5}$ | 1.2336e+3 $-$ | 5.5668e+5 $-$ | 1.6919e+6 $-$ | 5.1996e+3 $-$ | 5.8033e+6 $-$ | 1.5258e+5 $-$ | 7.2347e+2 $-$ | 1.2120e+2 $\approx$ | 2.0580e+2
(6.8051e+2) | (2.7471e+5) | (1.0907e+6) | (2.8736e+3) | (2.9157e+6) | (7.7532e+4) | (3.3376e+2) | (6.9918e+1) | (1.981e+2)
$f_{6}$ | 1.9375e+1 $-$ | 3.4820e+3 $-$ | 4.5314e+3 $-$ | 1.3163e+2 $-$ | 8.8923e+3 $-$ | 1.0266e+3 $-$ | 3.2250e+1 $-$ | 3.6250e+0 $-$ | 0.0000e+0
(3.9978e+0) | (1.1301e+3) | (8.0852e+2) | (3.2824e+1) | (3.1734e+3) | (4.3859e+2) | (7.7965e+0) | (2.1998e+0) | (0.000e+0)
$f_{7}$ | 9.1447e-2 $-$ | 8.2203e-1 $-$ | 1.1222e+0 $-$ | 8.2239e-2 $-$ | 5.2040e+1 $-$ | 5.4577e-2 $-$ | 5.8113e-2 $-$ | 8.4790e-2 $-$ | 2.0943e-2
(3.8031e-2) | (2.2065e-1) | (3.6673e-1) | (2.6534e-2) | (2.7550e+1) | (2.1989e-2) | (1.4618e-2) | (2.7932e-2) | (8.192e-3)
$f_{8}$ | -1.4025e+4 $\approx$ | -8.2308e+3 $-$ | -1.1684e+4 $-$ | -1.2799e+4 $-$ | -1.0531e+4 $-$ | -1.2285e+4 $-$ | -1.4023e+4 $\approx$ | -1.4575e+4 $+$ | -1.3852e+4
(5.0347e+2) | (1.2979e+3) | (5.6921e+2) | (2.8077e+2) | (8.2223e+2) | (4.2013e+2) | (8.5233e+1) | (2.0746e+2) | (6.717e+2)
$f_{9}$ | 1.1504e+1 $+$ | 1.5864e+2 $-$ | 2.7333e+2 $-$ | 2.2788e+2 $-$ | 2.6914e+2 $-$ | 1.0901e+2 $-$ | 1.6771e+2 $-$ | 2.3842e+1 $-$ | 1.4546e+1
(3.3186e+0) | (1.7533e+1) | (1.4085e+1) | (2.3562e+1) | (2.7653e+1) | (1.7125e+1) | (1.4089e+1) | (5.6385e+0) | (2.163e+0)
$f_{10}$ | 1.7023e+0 $-$ | 1.2329e+1 $-$ | 1.3329e+1 $-$ | 4.0097e+0 $-$ | 1.4671e+1 $-$ | 7.6077e+0 $-$ | 2.8370e+0 $-$ | 1.2297e+0 $-$ | 2.1621e-1
(3.1651e-1) | (1.3322e+0) | (8.5004e-1) | (5.6144e-1) | (7.3884e-1) | (1.1400e+0) | (2.1827e-1) | (3.8807e-1) | (5.119e-2)
$f_{11}$ | 1.1563e+0 $-$ | 2.8566e+1 $-$ | 4.6051e+1 $-$ | 2.1092e+0 $-$ | 1.1710e+2 $-$ | 1.0340e+1 $-$ | 1.2644e+0 $-$ | 1.0158e+0 $-$ | 6.7928e-1
(4.0584e-2) | (3.2736e+0) | (1.3901e+1) | (2.8414e-1) | (3.2432e+1) | (3.2046e+0) | (7.1606e-2) | (3.6145e-2) | (1.999e-1)
$f_{12}$ | 2.8786e+1 $-$ | 1.0624e+3 $-$ | 7.1756e+4 $-$ | 4.0254e+1 $-$ | 4.9839e+6 $-$ | 1.3627e+2 $-$ | 3.1790e+1 $-$ | 7.2029e+1 $-$ | 1.5847e+0
(1.5748e+1) | (1.3982e+3) | (6.5424e+4) | (5.5667e+0) | (5.9492e+6) | (3.5766e+1) | (9.1497e+0) | (2.4378e+1) | (1.976e+0)
$f_{13}$ | 7.1895e+0 $-$ | 4.3602e+5 $-$ | 2.1898e+6 $-$ | 1.0939e+1 $-$ | 1.2093e+7 $-$ | 2.1808e+3 $-$ | 4.3094e+0 $-$ | 4.3502e+0 $-$ | 1.0075e-1
(7.2480e+0) | (5.1304e+5) | (1.7255e+6) | (2.7540e+0) | (1.0145e+7) | (5.2706e+3) | (1.5238e+0) | (2.3093e+0) | (4.216e-2)
$+/-/\approx$ | 1/11/1 | 0/13/0 | 0/13/0 | 0/13/0 | 0/13/0 | 0/12/1 | 0/11/2 | 2/10/1 |

  * ”†”

’$+$’, ’$-$’ and ’$\approx$’ indicate that the result is significantly better,
significantly worse, and statistically similar to that obtained by AutoV.

To study the influence of the number of parameter sets $k$, Fig. 5 plots the
performance of the metaheuristic with operator $h_{3}$ and $k=1,5,10,15,20$ on
the eight benchmark problems, where $k=10$ leads to better overall performance
than the other settings. On the one hand, a small value of $k$ provides a few
different search behaviors, and thus leads to a low performance limit. On the
other hand, although a large value of $k$ provides many different search
behaviors, it leads to a large number of parameters that are difficult to be
optimized. As a consequence, $k=10$ is a proper setting for finding high-
performance operators.

TABLE VIII: Minimum Objective Values Obtained by Nine Metaheuristics on 13
Rotated Small-Scale Benchmark Problems with 30 Decision Variables. Best
Results are Highlighted.

Rotated |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
small-scale | GA [6] | PSO [10] | DE [7] | CMA-ES [8] | FEP [9] | CSO [37] | SHADE [38] | IMODE [39] | AutoV
problems [9] |  |  |  |  |  |  |  |  |
$f_{1}$ | 1.7644e+1 $-$ | 3.1210e+3 $-$ | 4.3448e+3 $-$ | 1.2422e+2 $-$ | 9.2516e+3 $-$ | 1.2894e+3 $-$ | 2.9700e+1 $-$ | 2.0650e+0 $-$ | 7.1625e-1
(5.3457e+0) | (8.2974e+2) | (8.3965e+2) | (2.9570e+1) | (3.9378e+3) | (6.7999e+2) | (1.0900e+1) | (1.0838e+0) | (5.403e-1)
$f_{2}$ | 1.9315e+1 $-$ | 3.7244e+1 $-$ | 4.9671e+1 $-$ | 1.3854e+1 $-$ | 8.1606e+1 $-$ | 1.4780e+1 $-$ | 1.8746e+1 $-$ | 1.5749e+0 $\approx$ | 1.5246e+0
(1.0234e+1) | (1.8184e+1) | (8.5749e+0) | (6.8676e+0) | (9.7066e+0) | (3.0538e+0) | (2.6580e+0) | (4.2566e-1) | (5.190e-1)
$f_{3}$ | 9.9581e+3 $-$ | 5.3590e+3 $-$ | 5.3369e+3 $-$ | 1.3929e+4 $-$ | 2.2100e+4 $-$ | 1.5239e+3 $\approx$ | 1.7530e+3 $\approx$ | 8.5084e+2 $+$ | 2.1193e+3
(2.6366e+3) | (1.4792e+3) | (1.3071e+3) | (3.9351e+3) | (5.0673e+3) | (9.7521e+2) | (4.7904e+2) | (2.5578e+2) | (7.960e+2)
$f_{4}$ | 1.4962e+1 $-$ | 2.4756e+1 $-$ | 3.6628e+1 $-$ | 7.8044e+0 $-$ | 4.6193e+1 $-$ | 1.2863e+1 $-$ | 8.1105e+0 $-$ | 5.5831e+0 $-$ | 3.2227e+0
(5.9860e+0) | (4.8244e+0) | (4.7317e+0) | (9.4482e-1) | (2.4146e+0) | (1.8236e+0) | (8.4983e-1) | (1.2362e+0) | (1.742e+0)
$f_{5}$ | 2.8300e+4 $-$ | 8.9721e+5 $-$ | 1.7625e+6 $-$ | 9.3046e+3 $-$ | 5.2830e+6 $-$ | 1.1350e+5 $-$ | 2.8513e+3 $-$ | 2.9329e+2 $\approx$ | 2.0175e+2
(4.8151e+4) | (3.4766e+5) | (7.4240e+5) | (6.0657e+3) | (2.3632e+6) | (6.7006e+4) | (3.4473e+3) | (1.6130e+2) | (1.540e+2)
$f_{6}$ | 2.3375e+1 $-$ | 3.6026e+3 $-$ | 5.0828e+3 $-$ | 1.1925e+2 $-$ | 8.6534e+3 $-$ | 1.2864e+3 $-$ | 3.2250e+1 $-$ | 1.3750e+1 $-$ | 2.5000e+0
(3.7009e+0) | (6.8520e+2) | (1.1280e+3) | (2.9011e+1) | (2.1819e+3) | (5.8652e+2) | (6.5629e+0) | (7.1464e+0) | (2.204e+0)
$f_{7}$ | 8.1742e-2 $-$ | 8.6616e-1 $-$ | 1.1051e+0 $-$ | 7.1672e-2 $-$ | 4.8419e+1 $-$ | 1.1136e-1 $-$ | 6.3319e-2 $-$ | 4.9230e-2 $-$ | 1.7558e-2
(2.9263e-2) | (4.3408e-1) | (6.6114e-1) | (3.4086e-2) | (2.8474e+1) | (5.2390e-2) | (1.6838e-2) | (2.2971e-2) | (5.800e-3)
$f_{8}$ | -7.6960e+3 $\approx$ | -6.4301e+3 $-$ | -5.2233e+3 $-$ | -7.6329e+3 $\approx$ | -7.4567e+3 $\approx$ | -7.5311e+3 $\approx$ | -5.5305e+3 $-$ | -7.1206e+3 $\approx$ | -7.3751e+3
(4.2092e+2) | (7.3448e+2) | (4.2305e+2) | (5.6710e+2) | (6.5244e+2) | (2.9851e+2) | (2.9295e+2) | (2.0444e+2) | (2.951e+2)
$f_{9}$ | 7.0727e+1 $-$ | 1.4357e+2 $-$ | 2.6278e+2 $-$ | 2.2499e+2 $-$ | 3.0943e+2 $-$ | 1.0129e+2 $-$ | 2.1434e+2 $-$ | 8.0586e+1 $-$ | 4.1289e+1
(1.4655e+1) | (1.9674e+1) | (1.3668e+1) | (2.5540e+1) | (1.3695e+1) | (1.1932e+1) | (1.3079e+1) | (1.7108e+1) | (1.364e+1)
$f_{10}$ | 2.9250e+0 $-$ | 1.1749e+1 $-$ | 1.3413e+1 $-$ | 4.3282e+0 $-$ | 1.5657e+1 $-$ | 7.3892e+0 $-$ | 3.0076e+0 $-$ | 2.3761e+0 $-$ | 9.6926e-1
(1.0910e-1) | (5.8613e-1) | (8.0644e-1) | (4.3019e-1) | (1.5809e+0) | (9.1288e-1) | (3.2128e-1) | (5.8902e-1) | (6.627e-1)
$f_{11}$ | 1.1551e+0 $-$ | 2.9789e+1 $-$ | 4.5249e+1 $-$ | 1.9904e+0 $-$ | 1.1321e+2 $-$ | 1.2782e+1 $-$ | 1.1978e+0 $-$ | 9.4447e-1 $-$ | 7.6236e-1
(2.4035e-2) | (8.3839e+0) | (8.3656e+0) | (2.2665e-1) | (2.2152e+1) | (4.1098e+0) | (7.2707e-2) | (7.3867e-2) | (1.343e-1)
$f_{12}$ | 3.5756e+1 $-$ | 1.5007e+3 $-$ | 1.9844e+5 $-$ | 3.7031e+1 $-$ | 1.2756e+6 $-$ | 1.3627e+2 $-$ | 3.4071e+1 $-$ | 7.2069e+1 $-$ | 1.4346e+0
(1.7714e+1) | (1.6007e+3) | (1.8928e+5) | (6.1593e+0) | (1.1497e+6) | (6.7341e+1) | (5.6141e+0) | (1.6143e+1) | (1.133e+0)
$f_{13}$ | 4.1967e+0 $-$ | 6.7210e+5 $-$ | 3.8989e+6 $-$ | 1.0194e+1 $-$ | 9.2196e+6 $-$ | 7.5595e+3 $-$ | 5.1567e+0 $-$ | 6.4904e+0 $-$ | 1.4557e-1
(3.5806e+0) | (4.4019e+5) | (2.7548e+6) | (3.2998e+0) | (7.2385e+6) | (2.0845e+4) | (1.6990e+0) | (1.5676e+0) | (6.230e-2)
$+/-/\approx$ | 0/12/1 | 0/13/0 | 0/13/0 | 0/12/1 | 0/12/1 | 0/11/2 | 0/12/1 | 1/9/3 |

  * ”†”

’$+$’, ’$-$’ and ’$\approx$’ indicate that the result is significantly better,
significantly worse, and statistically similar to that obtained by AutoV.

Furthermore, the influence of the benchmark problem on fitness evaluation is
studied. Table VI lists the performance of the metaheuristic with operator
$h_{3}$ optimized on four benchmark problems, including the Schwefel’s
Function 2.22 with a unimodal landscape, the Schwefel’s Function 2.21 with a
flat landscape, and the Ackley’s Function and Rastrigin’s Function with
multimodal landscapes. It can be found that the four metaheuristics obtain the
best performance on the benchmark problem for fitness evaluation, while they
obtain quite similar performance on the other problems. In the following
experiments, the metaheuristic with operator $h_{3}$ and $k=10$ optimized on
the Rastrigin’s Function is used as a representative of AutoV to be compared
with existing metaheuristics on various problems. Details of the variation
operator found by AutoV are presented in the following, where the components
with $r_{i}\sim\mathcal{N}(0,0)$ are ignored for simplicity:

$\displaystyle h(x_{1},x_{2},l,u)=$ (81) $\displaystyle\begin{aligned}
\left\\{\begin{aligned}
&(1-r_{2})x_{1}+r_{2}x_{2},&r_{2}\sim\mathcal{N}(0.4753,0.0103)\\\ &&{\rm if}\
p<0.219\\\ &(1-r_{2})x_{1}+r_{2}x_{2},&r_{2}\sim\mathcal{N}(0.0034,0.0006)\\\
&&{\rm if}\ 0.219\leq p<0.436\\\
&(1-r_{2})x_{1}+r_{2}x_{2},&r_{2}\sim\mathcal{N}(0.9999,0.0072)\\\ &&{\rm if}\
0.436\leq p<0.644\\\
&(1-r_{2})x_{1}+r_{2}x_{2},&r_{2}\sim\mathcal{N}(0.9988,0.0167)\\\ &&{\rm if}\
0.644\leq p<0.802\\\
&(1-r_{2})x_{1}+r_{2}x_{2},&r_{2}\sim\mathcal{N}(-0.0070,0.0056)\\\ &&{\rm
if}\ 0.802\leq p<0.960\\\
&(1-r_{2}-r_{3}-r_{4})x_{1}+r_{2}x_{2}+r_{3}l+r_{4}u,&r_{2}\sim\mathcal{N}(0.7293,0.8055)\\\
&&r_{3}\sim\mathcal{N}(0,0.0014)\\\ &&r_{4}\sim\mathcal{N}(0,0.0014)\\\ &&{\rm
if}\ 0.960\leq p<0.992\\\
&(1-r_{2}-r_{3}-r_{4})x_{1}+r_{2}x_{2}+r_{3}l+r_{4}u,&r_{2}\sim\mathcal{N}(0.0392,0.3185)\\\
&&r_{3}\sim\mathcal{N}(0,0.0351)\\\ &&r_{4}\sim\mathcal{N}(0,0.0351)\\\ &&{\rm
if}\ 0.992\leq p<0.997\\\
&(1-r_{2}-r_{3}-r_{4})x_{1}+r_{2}x_{2}+r_{3}l+r_{4}u,&r_{2}\sim\mathcal{N}(-0.9998,0.2478)\\\
&&r_{3}\sim\mathcal{N}(0,0.0011)\\\ &&r_{4}\sim\mathcal{N}(0,0.0011)\\\ &&{\rm
if}\ 0.997\leq p<0.998\\\
&(1-r_{2}-r_{3}-r_{4})x_{1}+r_{2}x_{2}+r_{3}l+r_{4}u,&r_{2}\sim\mathcal{N}(-0.8547,0.1174)\\\
&&r_{3}\sim\mathcal{N}(0,0.0452)\\\ &&r_{4}\sim\mathcal{N}(0,0.0452)\\\ &&{\rm
if}\ 0.998\leq p<0.999\\\
&(1-r_{2}-r_{3}-r_{4})x_{1}+r_{2}x_{2}+r_{3}l+r_{4}u,&r_{2}\sim\mathcal{N}(-0.5621,0.1059)\\\
&&r_{3}\sim\mathcal{N}(0,0.0003)\\\ &&r_{4}\sim\mathcal{N}(0,0.0003)\\\ &&{\rm
if}\ p\geq 0.999\\\ \end{aligned}\right.,\end{aligned}$

where $p$ is a uniformly distributed random value sampled in $[0,1]$. It can
be found from the definition that the operator contains ten functions whose
parameters obey different distributions, which provide a variety of search
behaviors that may be effective for different problems. It is worth noting
that the ten functions are selected with different probabilities, where the
operator uses only $x_{1}$ and $x_{2}$ in most cases but rarely uses the lower
bound $l$ and the upper bound $u$. This is similar to the genetic algorithm,
which performs the crossover between two parents with a large probability for
global search and performs the mutation on a single solution with a small
probability for local search, where the lower and upper bounds are involved in
the mutation operator (e.g., polynomial mutation [36]) but ignored in the
crossover operator (e.g., SBX [30]).

###  V-B Comparison on Small-Scale Benchmark Problems

Then, the proposed AutoV (i.e., the metaheuristic with operator $h_{3}$) is
compared with eight existing metaheuristics on 13 small-scale benchmark
problems with 30 decision variables, where the definitions of these problems
can be found in [9]. For all the compared metaheuristics, the population size
is set to 100 and the number of function evaluations is set to 10000. Table
VII shows the means and standard deviations of the minimum objective values
found by the nine metaheuristics, averaged over 30 runs. It can be found that
the proposed AutoV obtains the best overall performance, which gains the best
results on 9 out of 13 problems. The Wilcoxon rank sum test [40] with a
significance level of 0.05 is adopted to perform statistical analysis, where
AutoV is significantly better than GA, PSO, DE, CMA-ES, FEP, CSO, SHADE, IMODE
on 11, 13, 13, 13, 13, 12, 11, and 10 problems, respectively. Furthermore,
Fig. 6 depicts the convergence trajectories of the compared metaheuristics on
the Step Function and the Penalized Function, where AutoV converges faster
than the other metaheuristics.

(a)  Figure 6: Convergence trajectories of nine metaheuristics on two
problems. TABLE IX: Minimum Objective Values Obtained by Nine Algorithms on
the CEC’2013 Large-Scale Benchmark Problems with About 1000 Decision
Variables. Best Results are Highlighted.

CEC’2013 |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
large-scale | GA [6] | PSO [10] | DE [7] | CMA-ES [8] | FEP [9] | CSO [37] | SHADE [38] | IMODE [39] | AutoV
problems [41] |  |  |  |  |  |  |  |  |
$f_{1}$ | 1.1824e+9 $-$ | 1.7274e+11 $-$ | 1.1820e+11 $-$ | 1.3230e+10 $-$ | 6.2225e+10 $-$ | 1.6506e+11 $-$ | 9.2604e+8 $-$ | 1.5672e+10 $-$ | 7.5437e+8
(1.5088e+8) | (8.2047e+9) | (5.7453e+9) | (5.9476e+8) | (8.8355e+9) | (6.4000e+9) | (9.2525e+7) | (1.8277e+9) | (3.599e+7)
$f_{2}$ | 7.0854e+3 $+$ | 4.8900e+4 $-$ | 4.2656e+4 $-$ | 1.1326e+4 $+$ | 8.4194e+4 $-$ | 4.4704e+4 $-$ | 1.8598e+4 $-$ | 2.4830e+4 $-$ | 1.2736e+4
(5.5474e+2) | (1.6515e+3) | (9.8044e+2) | (2.6942e+2) | (5.6859e+2) | (7.5103e+2) | (1.1309e+3) | (6.9359e+2) | (3.996e+2)
$f_{3}$ | 1.9790e+1 $-$ | 2.1096e+1 $-$ | 2.1011e+1 $-$ | 2.1334e+1 $-$ | 2.1457e+1 $-$ | 2.0977e+1 $-$ | 1.7116e+1 $-$ | 2.0013e+1 $-$ | 7.6972e+0
(1.0792e-1) | (3.6567e-2) | (6.4041e-3) | (2.0615e-2) | (8.6678e-3) | (1.9539e-2) | (2.1804e-1) | (1.6968e-2) | (2.709e-1)
$f_{4}$ | 7.0214e+11 $-$ | 3.8264e+12 $-$ | 2.8759e+11 $-$ | 1.8046e+12 $-$ | 7.4290e+11 $-$ | 2.0687e+12 $-$ | 4.1680e+10 $+$ | 3.5605e+11 $-$ | 1.7434e+11
(3.6361e+11) | (1.0018e+12) | (4.4452e+10) | (2.2407e+11) | (4.1395e+11) | (5.9504e+11) | (7.7120e+9) | (3.1439e+10) | (3.356e+10)
$f_{5}$ | 8.0929e+6 $-$ | 3.3051e+7 $-$ | 1.6648e+7 $-$ | 1.1814e+7 $-$ | 2.0252e+7 $-$ | 2.7747e+7 $-$ | 6.6297e+6 $-$ | 1.0851e+7 $-$ | 5.5157e+6
(1.1105e+6) | (3.3554e+6) | (1.2160e+6) | (3.5863e+5) | (1.9684e+6) | (2.3903e+6) | (8.5226e+5) | (6.3028e+5) | (7.817e+5)
$f_{6}$ | 7.6705e+5 $-$ | 1.0156e+6 $-$ | 9.6778e+5 $-$ | 1.0725e+6 $-$ | 7.7643e+5 $-$ | 9.9765e+5 $-$ | 1.0802e+5 $-$ | 2.7032e+5 $-$ | 8.3894e+4
(4.2587e+4) | (6.2467e+3) | (1.9218e+4) | (1.4402e+3) | (6.3779e+4) | (5.1452e+3) | (1.1582e+4) | (4.2843e+4) | (9.288e+3)
$f_{7}$ | 7.8310e+9 $-$ | 4.9258e+13 $-$ | 2.3870e+12 $-$ | 1.1084e+10 $-$ | 2.8185e+11 $-$ | 2.2211e+13 $-$ | 8.9123e+8 $\approx$ | 1.3523e+10 $-$ | 8.4166e+8
(2.9908e+9) | (2.3576e+13) | (7.2173e+11) | (9.4838e+9) | (1.4497e+11) | (8.4947e+12) | (1.9320e+8) | (1.6130e+9) | (7.997e+7)
$f_{8}$ | 2.9461e+16 $-$ | 1.6998e+17 $-$ | 1.2728e+13 $+$ | 4.7204e+16 $-$ | 6.1725e+15 $\approx$ | 3.2323e+16 $-$ | 3.7727e+13 $+$ | 5.4599e+15 $\approx$ | 4.4372e+15
(1.4136e+16) | (9.1594e+16) | (1.9676e+12) | (6.0409e+15) | (3.8102e+15) | (3.0414e+16) | (3.0187e+13) | (2.5604e+15) | (1.770e+15)
$f_{9}$ | 7.5728e+8 $-$ | 2.5298e+9 $-$ | 1.2725e+9 $-$ | 8.4502e+8 $-$ | 1.6577e+9 $-$ | 2.0244e+9 $-$ | 5.8016e+8 $-$ | 7.7668e+8 $-$ | 4.4601e+8
(1.2373e+8) | (2.0469e+8) | (6.3827e+7) | (4.5495e+7) | (1.3394e+8) | (2.1854e+8) | (1.1320e+8) | (4.1594e+7) | (8.503e+7)
$f_{10}$ | 4.4221e+7 $-$ | 8.7041e+7 $-$ | 2.1046e+6 $-$ | 9.6156e+7 $-$ | 1.7440e+7 $-$ | 7.6458e+7 $-$ | 1.2506e+6 $-$ | 7.8125e+5 $\approx$ | 4.8728e+5
(1.3077e+7) | (2.6524e+6) | (6.4391e+5) | (4.6215e+5) | (3.3440e+6) | (5.2418e+6) | (1.4290e+4) | (2.8894e+5) | (4.870e+5)
$f_{11}$ | 7.4972e+11 $-$ | 3.0398e+15 $-$ | 1.1061e+14 $-$ | 4.0919e+11 $-$ | 1.1334e+13 $-$ | 8.9529e+14 $-$ | 8.7764e+9 $+$ | 4.4608e+11 $-$ | 1.6943e+11
(4.5447e+11) | (1.4680e+15) | (9.8496e+13) | (1.0871e+11) | (5.1094e+12) | (3.6702e+14) | (4.5005e+9) | (1.6708e+11) | (1.144e+11)
$f_{12}$ | 2.0599e+10 $-$ | 1.9723e+12 $-$ | 1.5925e+12 $-$ | 1.3548e+10 $-$ | 2.6665e+12 $-$ | 1.6723e+12 $-$ | 3.3481e+10 $-$ | 6.4235e+11 $-$ | 7.6900e+8
(3.9186e+9) | (8.1522e+10) | (3.3820e+10) | (1.3457e+10) | (1.0994e+11) | (2.5065e+10) | (5.8270e+9) | (1.9268e+10) | (2.144e+7)
$f_{13}$ | 4.0277e+10 $-$ | 1.9190e+15 $-$ | 5.4725e+13 $-$ | 8.3330e+10 $-$ | 3.0193e+12 $-$ | 1.0082e+15 $-$ | 1.0448e+10 $+$ | 6.6059e+10 $-$ | 1.7051e+10
(1.4024e+10) | (8.0609e+14) | (1.7016e+13) | (1.4211e+10) | (1.6034e+12) | (6.2874e+14) | (2.3651e+9) | (9.5957e+9) | (4.590e+9)
$f_{14}$ | 1.0031e+12 $-$ | 4.5075e+15 $-$ | 1.2185e+14 $-$ | 9.6175e+11 $-$ | 1.2889e+13 $-$ | 1.1075e+15 $-$ | 7.7424e+10 $+$ | 7.2323e+11 $-$ | 2.3734e+11
(3.8285e+11) | (3.3588e+15) | (7.0856e+13) | (4.0377e+11) | (4.1028e+12) | (4.3764e+14) | (1.3029e+10) | (8.4437e+10) | (7.915e+10)
$f_{15}$ | 5.7586e+7 $-$ | 1.8722e+15 $-$ | 3.4395e+14 $-$ | 9.9936e+8 $-$ | 8.8170e+13 $-$ | 8.7705e+14 $-$ | 3.0449e+7 $\approx$ | 6.1844e+12 $-$ | 2.7840e+7
(1.5169e+7) | (2.4943e+14) | (4.9325e+13) | (1.1950e+9) | (2.5655e+13) | (1.2736e+14) | (4.6006e+6) | (1.4191e+12) | (2.630e+6)
$+/-/\approx$ | 1/14/0 | 0/15/0 | 1/14/0 | 1/14/0 | 0/14/1 | 0/15/0 | 5/8/2 | 0/13/2 |

  * ”†”

’$+$’, ’$-$’ and ’$\approx$’ indicate that the result is significantly better,
significantly worse, and statistically similar to that obtained by AutoV.

While the 13 benchmark problems are already translated and scaled, they are
further rotated by randomly generated orthogonal matrices and challenge the
nine metaheuristics. As can be seen from the experimental results listed in
Table VIII, the superiority of the proposed AutoV becomes more significant,
where AutoV outperforms the other metaheuristics on 11 out of 13 problems. As
a consequence, the superiority of AutoV over some classical and state-of-the-
art metaheuristics can be verified. Besides, it also implies that AutoV is
more effective than the approaches based on the recommendation and combination
of existing metaheuristics, whose performance can hardly go beyond the best
existing metaheuristic on each problem.

###  V-C Comparison on Large-Scale Benchmark Problems

Lastly, the proposed AutoV and the eight existing metaheuristics are compared
on the 15 CEC’2013 large-scale benchmark problems [41]. These benchmark
problems contain approximately 1000 decision variables and a variety of
landscape functions, transformations, and interactions between variables,
posing stiff challenges to general metaheuristics. For all the compared
algorithms, the population size is set to 100 and the number of function
evaluations is set to 120000.

Table IX presents the means and standard deviations of the minimum objective
values found by the compared metaheuristics, averaged over 30 runs. It can be
observed from the statistical results that the proposed AutoV also exhibits
better overall performance than the others on the large-scale benchmark
problems, achieving the best results on 9 out of 15 problems. It is noteworthy
that although SHADE and IMODE suggest many complex search strategies for the
combination of multiple operators and adaptation of parameters, they are still
underperformed by AutoV that only contains a simple operator designed
automatically. Therefore, the proposed AutoV offers bright prospects to the
design of metaheuristics, which can potentially replace the laborious manual
design process.

##  VI Conclusions

To reduce the human expertise in designing metaheuristics, this paper has
analyzed the importance of translation, scale, and rotation invariance to the
robustness of operators, and deduced the generic form of translation, scale,
rotation invariant operators. Based on the deduced generic form, this paper
has proposed a principled approach to search for high-performance operators
automatically. In contrast to the automated design approaches based on
existing metaheuristics, the proposed approach does not rely on any existing
techniques, and can obtain competitive performance over some state-of-the-art
metaheuristics on complex and large-scale optimization problems.

The experimental results have demonstrated the effectiveness and potential of
the automated design of variation operators, and further research on this
topic is highly desirable. Firstly, it is reasonable to search for high-
performance operators based on more complex functions (e.g., include more
parents and consider the update of velocity), where more effective operators
are expected to be found. Secondly, since the proposed approach searches for
operators according to their performance on a benchmark problem, it is
reasonable to adopt more representative and practical problems to find high-
performance operators for general problems or specific types of problems.
Thirdly, it is interesting to develop novel approaches for automatically
designing selection strategies for metaheuristics.

## References

  * [1] M. R. Alam, K. S. Lee, M. Rahman, and Y. F. Zhang, “Process planning optimization for the manufacture of injection moulds using a genetic algorithm,” _International Journal of Computer Integrated Manufacturing_ , vol. 16, no. 3, pp. 181–191, 2003.
  * [2] J. Y. Potvin and S. Bengio, “The vehicle routing problem with time windows part II: Genetic search,” _INFORMS Journal on Computing_ , vol. 8, no. 2, pp. 165–172, 1996.
  * [3] Y. Tian, X. Su, Y. Su, and X. Zhang, “EMODMI: A multi-objective optimization based method to identify disease modules,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , 2020.
  * [4] K. J. Oh, T. Y. Kim, and S. Min, “Using genetic algorithm to support portfolio optimization for index fund management,” _Expert Systems with Applications_ , vol. 28, pp. 371–379, 2005.
  * [5] H. Robbins and S. Monro, “A stochastic approximation method,” _The Annals of Mathematical Statistics_ , pp. 400–407, 1951.
  * [6] J. H. Holland, _Adaptation in Natural and Artificial Systems_. MIT Press, 1992.
  * [7] R. Storn and K. Price, “Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces,” _Journal of Global Optimization_ , vol. 11, no. 4, pp. 341–359, 1997.
  * [8] N. Hansen and A. Ostermeier, “Completely derandomized self-adaptation in evolution strategies,” _Evolutionary Computation_ , vol. 9, no. 2, pp. 159–195, 2001.
  * [9] X. Yao, Y. Liu, and G. Lin, “Evolutionary programming made faster,” _IEEE Transactions on Evolutionary Computation_ , vol. 3, no. 2, pp. 82–102, 1999.
  * [10] R. Eberhart and J. Kennedy, “A new optimizer using particle swarm theory,” in _Proceedings of the 6th International Symposium on Micro Machine and Human Science_ , 1995, pp. 39–43.
  * [11] K. Socha and M. Dorigo, “Ant colony optimization for continuous domains,” _European Journal of Operational Research_ , vol. 185, no. 3, pp. 1155–1173, 2008.
  * [12] D. Karaboga, “An idea based on honey bee swarm for numerical optimization,” Erciyes University, Tech. Rep. tr06, 2005.
  * [13] J. Galletly, _Evolutionary algorithms in theory and practice: Evolution strategies, evolutionary programming, genetic algorithms_. Kybernetes, 1998.
  * [14] T. Okabe, Y. Jin, and B. Sendhoff, “Theoretical comparisons of search dynamics of genetic algorithms and evolution strategies,” in _Proceedings of the 2005 IEEE Congress on Evolutionary Computation_ , 2005, pp. 382–389.
  * [15] S. Li, Y. Li, and Y. Lin, _Intelligent optimization algorithms and emergent computation_. Beijing: Tsinghua University Press, 2019.
  * [16] A. E. Eiben and C. A. Schippers, “On evolutionary exploration and exploitation,” _Fundamenta Informaticae_ , vol. 35, pp. 1–16, 1998.
  * [17] Y. Su, N. Guo, Y. Tian, and X. Zhang, “A non-revisiting genetic algorithm based on a novel binary space partition tree,” _Information Sciences_ , vol. 512, pp. 661–674, 2020.
  * [18] H. Li and Q. Zhang, “Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II,” _IEEE Transactions on Evolutionary Computation_ , vol. 13, no. 2, pp. 284–302, 2009.
  * [19] T. Rodemann, “Industrial portfolio management for many-objective optimization algorithms,” in _Proceedings of the 2018 IEEE Congress on Evolutionary Computation_ , 2018.
  * [20] C. A. Coello Coello and M. S. Lechuga, “MOPSO: A proposal for multiple objective particle swarm optimization,” in _Proceedings of the 2002 IEEE Congress on Evolutionary Computation_ , vol. 2, 2002, pp. 1051–1056.
  * [21] M. Jebalia, A. Auger, and N. Hansen, “Log-linear convergence and divergence of the scale-invariant (1+1)-ES in noisy environments,” _Algorithmica_ , vol. 59, no. 3, pp. 425–460, 2011.
  * [22] F. Caraffini and F. Neri, “A study on rotation invariance in differential evolution,” _Swarm and Evolutionary Computation_ , vol. 50, p. 100436, 2019\.
  * [23] G. Karafotias, M. Hoogendoorn, and A. E. Eiben, “Parameter control in evolutionary algorithms: Trends and challenges,” _IEEE Transactions on Evolutionary Computation_ , vol. 19, no. 2, pp. 167–187, 2015.
  * [24] C. Huang, Y. Li, and X. Yao, “A survey of automatic parameter tuning methods for metaheuristics,” _IEEE Transactions on Evolutionary Computation_ , 2019, in press.
  * [25] Y. Tian, S. Peng, X. Zhang, T. Rodemann, K. C. Tan, and Y. Jin, “A recommender system for metaheuristic algorithms for continuous optimization based on deep recurrent neural networks,” _IEEE Transactions on Artificial Intelligence_ , vol. 1, no. 1, pp. 5–18, 2020.
  * [26] L. C. Bezerra, M. López-Ibánez, and T. Stützle, “Automatic component-wise design of multiobjective evolutionary algorithms,” _IEEE Transactions on Evolutionary Computation_ , vol. 20, no. 3, pp. 403–417, 2016\.
  * [27] L. P. M. Birattari, T. Stützle and K. Varrentrapp, “A racing algorithm for configuring metaheuristics,” in _Proceedings of the 2002 Conference on Genetic and Evolutionary Computation_ , 2002, pp. 11–18.
  * [28] F. Rosenblatt, “The perceptron: a probabilistic model for information storage and organization in the brain,” _Psychological Review_ , vol. 65, no. 6, pp. 386–408, 1958.
  * [29] A. Fialho, M. Schoenauer, and M. Sebag, “Toward comparison-based adaptive operator selection,” in _Proceedings of the 2010 Conference on Genetic and Evolutionary Computation_ , 2002, pp. 767–774.
  * [30] K. Deb and R. B. Agrawal, “Simulated binary crossover for continuous search space,” _Complex Systems_ , vol. 9, no. 4, pp. 115–148, 1995.
  * [31] C. Igel, N. Hansen, and S. Roth, “Covariance matrix adaptation for multi-objective optimization,” _Evolutionary computation_ , vol. 15, no. 1, pp. 1–28, 2007.
  * [32] N. Hansen, R. Ros, N. Mauny, M. Schoenauer, and A. Auger, “Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems,” _Applied Soft Computing_ , vol. 11, no. 8, pp. 5755–5769, 2011.
  * [33] K. A. Hoffmann and S. T. Chiang, _Computational Fluid Dynamics Volume I_. Engineering Education System, 2000\.
  * [34] L. Pan, W. Xu, L. Li, C. He, and R. Cheng, “Adaptive simulated binary crossover for rotated multi-objective optimization,” _Swarm and Evolutionary Computation_ , vol. 60, p. 100759, 2021.
  * [35] D. B. Fogel, “Asymptotic convergence properties of genetic algorithms and evolutionary programming: Analysis and experiments,” _Cybernetics and Systems_ , vol. 25, no. 3, pp. 389–407, 1994.
  * [36] K. Deb and M. Goyal, “A combined genetic adaptive search (GeneAS) for engineering design,” _Computer Science and Informatics_ , vol. 26, no. 4, pp. 30–45, 1996.
  * [37] R. Cheng and Y. Jin, “A competitive swarm optimizer for large scale optimization,” _IEEE Transactions on Cybernetics_ , vol. 45, no. 2, pp. 191–204, 2014.
  * [38] R. Tanabe and A. Fukunaga, “Success-history based parameter adaptation for differential evolution,” in _Proceedings of the IEEE Congress on Evolutionary Computation_ , 2013.
  * [39] K. M. Sallam, S. M. Elsayed, R. K. Chakrabortty, and M. J. Ryan, “Improved multi-operator differential evolution algorithm for solving unconstrained problems,” in _Proceedings of the IEEE Congress on Evolutionary Computation_ , 2020.
  * [40] J. Derrac, S. Garcia, D. Molina, and F. Herrera, “A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms,” _Swarm and Evolutionary Computation_ , vol. 1, no. 1, pp. 3–18, 2011.
  * [41] X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, “Benchmark functions for the CEC’2013 special session and competition on large-scale global optimization,” _Gene_ , vol. 7, no. 33, 2013.