# Exploration Enhancement of Nature-Inspired Swarm-based Optimization
Algorithms

Kwok Pui Choi,Enzio Hai Hong Kam,Tze Leung Lai, Xin T. Tong, and Weng Kee Wong
K. Choi is with Department of Statistics and Applied Probability, National
University of SingaporeE. Kam is with the School of Computing, National
University of Singapore T. L. Lai is with Department of Statistics, Stanford
UniversityX. T. Tong is with Department of Mathematics, National University of
Singapore W. K. Wong is with Department of Biostatistics, University of
California, Los Angeles Manuscript received ???, ????; revised ???, ????. This
work is supported in part by Singapore MOE Academic Research Fund
R-155-000-222-114, .

###### Abstract

Nature-inspired swarm-based algorithms have been widely applied to tackle
high-dimensional and complex optimization problems across many disciplines.
They are general purpose optimization algorithms, easy to use and implement,
flexible and assumption-free. A common drawback of these algorithms is
premature convergence and the solution found is not a global optimum. We
provide sufficient conditions for an algorithm to converge almost surely
(a.s.) to a global optimum. We then propose a general, simple and effective
strategy, called Perturbation-Projection (PP), to enhance an algorithm’s
exploration capability so that our convergence conditions are guaranteed to
hold. We illustrate this approach using three widely used nature-inspired
swarm-based optimization algorithms: particle swarm optimization (PSO), bat
algorithm (BAT) and competitive swarm optimizer (CSO). Extensive numerical
experiments show that each of the three algorithms with the enhanced PP
strategy outperforms the original version in a number of notable ways.

###### Index Terms:

Bat algorithm, Competitive swarm optimizer, Convergence to the global optimum,
Exploration, Exploitation, Particle swarm optimization.

##  I Introduction

Over the past decades, Swarm Intelligence has and continues to inspire a
steadily rising numbers of nature-inspired swarm-based algorithms in
optimizing high-dimensional complex cost functions, including those that do
not have analytic forms. Examples of swarm-based algorithms include particle
swarm optimization (PSO), competitive swarm optimizer (CSO), ant colony
optimization, bat algorithm (BAT) etc. Swarm-based optimization algorithms are
mainly motivated by nature or animal behavior and then thoughtfully formulated
into an algorithm that iterates to the optimum according to a couple of
equations. Generally, these algorithms are easy to code and implement, and do
not need gradient information or other sophisticated mathematical techniques
to work well. Computer codes are widely and freely available, which have
undoubtedly helped fuel numerous and various applications of these algorithms
to tackle many different types of complex real-world optimization problems.
Documentation of their effectiveness is widespread resulting in their meteoric
applications in an increasing number of sectors in industry and in academia
[1, 2].

There are known shortcomings of nature-inspired swarm-based optimization
algorithms. We highlight two: (1) They require effective tuning of their
parameters to achieve optimal performance [3] and (2) they often suffer from
premature convergence to a local optimum of the cost function, especially when
the cost function is high-dimensional and multi-modal. Recently, [4] proposed
a general strategy to tackle the first shortcoming using statistical tools.

In this paper, we propose an innovative and simple strategy to address the
latter problem by modifying a nature-inspired swarm-based algorithm to ensure
its almost sure (a.s) convergence to a global optimum. The proposed
methodology is general, as long as certain unrestrictive technical conditions
are met, and we demonstrate the techniques using three popular swarm-based
algorithms mentioned above. Theorem II.1 in Section II shows that if a swarm-
based algorithm satisfy conditions (C1)–(C3), the algorithm is guaranteed to
converge a.s. to a global solution. We also provide comments on possible
extension of the proposed methodology to evolutionary algorithms, like
differential evolution.

We next propose a general, yet simple and effective strategy, Perturbation-
Projection (PP), to modify an algorithm, if necessary, to ensure conditions
(C1) and (C2) to hold. In Section III, we demonstrate how the PP strategy is
applied to the original PSO, BAT and CSO algorithms with appropriate
modifications to meet the required conditions for almost sure convergence to
the global optimum. We denote these modified algorithms of their counterparts
by mPSO, mBAT and mCSO. When it is necessary to further improve the
exploitation ability of an algorithm but not at the expense of its exploration
ability, we introduce a Heterogeneous Perturbation-Projection (HPP) strategy
in Section IV to satisfy conditions (C1)–(C2). In Section V, we conduct
extensive numerical experiments using many commonly used test functions and
show that the modified algorithms remain very competitive with the original
algorithms. We also conduct tuning parameter analysis for the PP and HPP
strategies. The paper ends with conclusion and related future work.

In the literature, there are many different modifications of the original
algorithms or new techniques to address the premature convergence issue for
swarm-based algorithms; see, for example [5, 6, 7, 8]. These methods tend to
either apply to only one specific algorithm, e.g. PSO, instead of a class of
algorithms, or, require additional detection or tuning mechanisms that can be
complicated to implement. A commonality among these works is that they do not
have rigorous mathematical theory to support the modifications for improved
convergence to the global optimum. The distinctive features of our proposed
modifications are that they are supported by mathematical theory, simple to
implement and applicable to any swarm-based algorithm.

##  II Stochastic enhancement of an algorithm’s exploration

In this section, we introduce a simple and effective strategy to modify a
nature-inspired swarm-based algorithm to ensure the modified algorithm
converges almost surely (a.s.) to a global optimum. The cost function can be a
high-dimensional, non-differentiable, non-separable or non-convex with
multiple local or global optima. The key idea behind the proposed algorithms
is simple; we incorporate an additional stochastic component with an
appropriate noise level to a nature-inspired swarm-based algorithm to enable
it to escape from a local minimum. While this idea may not be entirely new,
the distinguishing feature of proposed modification is its generality and the
modified algorithm is guaranteed to converge a.s. to a global optimum.

Without loss of generality, we assume that we have a minimization problem as
$\min_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})=-\max_{\mathbf{x}\in\mathcal{D}}\\{-f(\mathbf{x})\\}$.
Here $f(\mathbf{x})$ is a user-specified real-valued cost function defined on
a given compact subset $\mathcal{D}$ of $\mathbb{R}^{d}$. For many real-world
applications, there are physical and budgetary constraints and we assume that
these constraints have been appropriately translated into the compact search
space. For many other problems, like those in machine learning, regularization
terms, sometimes denoted by $\rho\|\mathbf{x}\|_{1}$ (or
$\rho\|\mathbf{x}\|^{2}_{2}$) are often added to the nonnegative loss/cost
function $f$ to prevent over-fitting, where the regularization parameter
$\rho$ can often be selected through various cross validation methods [9, 10].
In such situations, we minimize the function
$F(\mathbf{x})=f(\mathbf{x})+\rho\|\mathbf{x}\|_{1}$ (or
$f(\mathbf{x})+\rho\|\mathbf{x}\|^{2}_{2}$) over $\mathcal{D}$, which is the
closed $L_{1}$-ball (or $L_{2}$-ball) centered at the origin with radius
$f(\mathbf{0})/\rho$. This is because for $\mathbf{x}\notin\mathcal{D}$,
$F(\mathbf{x})>F(\mathbf{0})$ so the global minima are in $\mathcal{D}$.

To solve the optimization problem, we resort to nature-inspired swarm-based
algorithms, which are increasingly used to solve complex and high-dimensional
optimization problems. They generally employ an evolutionary-like or swarm-
based strategy to search for an optimum by first randomly generating a user-
specified set of candidate solutions. Depending on the particular swarm-based
algorithm, candidate solutions are called agents or particles. These
algorithms embrace common principles that include good exploration and
exploitation abilities and they have stochastic components and tuning
parameters. For example, for a swarm-based type of algorithm, it can be
typically described as follows.

At time or iteration $t=1,2,\ldots$, let $\mathbf{x}_{i}(t)$ represent the
$i$-th particle’s position or candidate solution of the optimization problem.
Let $\\{\mathbf{x}_{i}(t):1\leq i\leq n\\}$ denote the collection of candidate
solutions of the swarm of size $n$. In addition, there is auxiliary
information from each member of the swarm and we denote them collectively by
$\\{\mathbf{v}_{i}(t):1\leq i\leq n\\}$. For example, in PSO, the auxiliary
information of each particle is its velocity with which it flies to the
current position. We denote them by

$\mathbf{x}(t)=(\mathbf{x}_{1}(t),\ldots,\mathbf{x}_{n}(t))~{}\text{and}~{}\mathbf{v}(t)=(\mathbf{v}_{1}(t),\ldots,\mathbf{v}_{n}(t)).$

Suppose the algorithm has a stochastic update rule of the following form: for
agent $i$

$\displaystyle\mathbf{x}_{i}(t+1)$ $\displaystyle=$
$\displaystyle\Phi_{i}(\mathbf{x}(t),\mathbf{v}(t),f),$
$\displaystyle\mathbf{v}_{i}(t+1)$ $\displaystyle=$
$\displaystyle\Psi_{i}(\mathbf{x}(t),\mathbf{v}(t),f),$

for $1\leq i\leq n$, where $\Phi_{i},\Psi_{i}$ are some functions of which the
outputs can be random. At either a pre-selected deterministic time or a
stopping time $T$, the algorithm is terminated and the best solution
$\min_{1\leq i\leq n}f(\mathbf{x}_{i}(T))$ is output.

A key desirable property of a global optimization algorithm is its ability to
identify a global optimum eventually. To achieve convergence to the global
optimum, Törn and Žilinskas [11] remark that the algorithm has to be able to
explore densely in the search space. Inspired by their remark, we formulate
the following conditions on the algorithm’s agents,
$\mathbf{x}_{1}(t),\ldots,\mathbf{x}_{n}(t)$.

  * (C1)

At any iteration, all agents stay in $\mathcal{D}$, i.e.,
$\mathbf{x}_{i}(t)\in\mathcal{D}$ for $1\leq i\leq n$ and $t=0,1,\ldots,T$.

  * (C2)

Conditional on the information up to the current iteration $t$, there exists
an agent at the next iteration, $t+1$, that has a positive probability to
explore a region of improvement in $\mathcal{D}$. Specifically, there exists
an $\alpha>0$, independent of $t$, such that for any $d$-dimensional ball $B$
in $\mathcal{D}$ satisfying

$f(\mathbf{y})\leq\min_{1\leq i\leq n}f(\mathbf{x}_{i}(t)),\quad\mbox{ for all
}\mathbf{y}\in B,$

the following holds

$\max_{1\leq i\leq n}\mathbb{P}_{t}(\mathbf{x}_{i}(t+1)\in B)\geq\alpha|B|,$

with $|B|$ denoting the volume of $B$ and $\mathbb{P}_{t}$ the conditional
probability with information at the $t$-th iteration.

  * (C3)

The algorithm is time improving, i.e.,

$\min_{1\leq i\leq n}f(\mathbf{x}_{i}(t+1))\leq\min_{1\leq i\leq
n}f(\mathbf{x}_{i}(t)),\quad a.s.$

######  Theorem II.1.

Suppose $f$ is a continuous function defined on a given compact subset
$\mathcal{D}$ of $\mathbb{R}^{d}$ with at least one minimizer in the interior
of $\mathcal{D}$. Any swarm-based algorithm that satisfies (C1)–(C3) will
converge to the global minimum almost surely (a.s). In other words, if the
swarm-based algorithm has $n$ agents, $\lim_{t\to\infty}\min_{1\leq i\leq
n}f(\mathbf{x}_{i}(t))=\min_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})$ a.s..

###### Proof.

Let an interior point $\mathbf{x}_{*}\in\mathcal{D}$ denote a minimizer of
$f$. Without loss of generality, we assume $f(\mathbf{x}_{*})=0$. By (C3),
$m(t):=\min_{1\leq i\leq n}f(\mathbf{x}_{i}(t))$ is a non-increasing sequence,
thus it suffices to show that $m(t)\stackrel{{\scriptstyle P}}{{\to}}0$ as
$t\to\infty$. Given a fixed $\epsilon>0$, by continuity of $f$, there exists
$\delta>0$ such that $0\leq f(\mathbf{x})\leq\epsilon$ for $\mathbf{x}\in
B(\mathbf{x}_{*},\delta):=\\{\mathbf{x}:\|\mathbf{x}-\mathbf{x}_{*}\|\leq\delta\\}$.
Since $\mathbf{x}_{*}$ is an interior point,
$B(\mathbf{x}_{*},\delta)\subset\mathcal{D}$ for a sufficiently small
$\delta$. By (C2), at each iteration, there is an agent $\mathbf{x}_{i}$ such
that either $f(\mathbf{x}_{i}(t))\leq\epsilon$ or

$\mathbb{P}_{t}\left(\mathbf{x}_{i(t+1)}(t+1)\in
B(\mathbf{x}_{*},\delta)\right)\geq\alpha|B(\mathbf{x}_{*},\delta)|=:\alpha(\epsilon).$

If $f(\mathbf{x}_{i}(t))\leq\epsilon$, then $m(t)\leq\epsilon$. If
$\mathbf{x}_{i(t+1)}(t+1)\in B(\mathbf{x}^{*},\delta)$, then $m(t+1)\leq
f(\mathbf{x}_{i(t+1)}(t+1))\leq\epsilon$. By the tower property of conditional
expectation,

$\displaystyle\mathbb{P}(m(t)>\epsilon)$
$\displaystyle\leq\mathbb{P}(\mathbf{x}_{i(s+1)}(s+1)\notin
B(\mathbf{x}^{*},\delta),\text{ for }s=1,\ldots,t-1)$
$\displaystyle=\mathbb{E}\prod_{s=1}^{t-1}\mathbb{P}_{s}(\mathbf{x}_{i(s+1)}(s+1)\notin
B(\mathbf{x}^{*},\delta))$
$\displaystyle<\mathbb{E}\prod_{s=1}^{t-1}(1-\alpha(\epsilon))=(1-\alpha(\epsilon))^{t-1}.$

Therefore, $m(t)\stackrel{{\scriptstyle P}}{{\to}}0$ as $t\to\infty$. ∎

Theorem II.1 shows that any swarm-based algorithm satisfying (C1)–(C3) is
guaranteed to converge to a global optimal solution eventually. Many swarm-
based algorithms satisfy (C3), but not all swarm-based algorithms stay within
the search space $\mathcal{D}$ in their iterations thus violating condition
(C1) or they are exploratory in the sense of condition (C2). However, Section
II-A below shows that, with slight modifications of an algorithm, we can
ensure the modified algorithm satisfies (C1)–(C2).

###  II-A Perturbation-projection strategy

A first step is to ensure all agents stay in $\mathcal{D}$. For this purpose,
we introduce the projection onto $\mathcal{D}$ ensuring the agents stay inside
the search space, $\mathcal{D}$, a compact subset in $\mathbb{R}^{d}$:

$\chi_{\mathcal{D}}(\mathbf{x})=\arg\min_{\mathbf{y}\in\mathcal{D}}\|\mathbf{y}-\mathbf{x}\|_{2},$

where $\|\mathbf{y}-\mathbf{x}\|_{2}$ denotes the Euclidean distance between
$\mathbf{x}$ and $\mathbf{y}$. In other words,
$\chi_{\mathcal{D}}(\mathbf{x})$ is a point in $\mathcal{D}$ that is closest
to $\mathbf{x}$. When ties occur, they can be broken arbitrarily. By
definition, $\chi_{\mathcal{D}}(\mathbf{x})=\mathbf{x}$ if
$\mathbf{x}\in\mathcal{D}$. In general, adding the projection to an
algorithm’s iterations reduces the agents’ distances to the optimal solution
[10]. In particular, if $\mathbf{x}_{*}\in\mathcal{D}$ and $\mathcal{D}$ is
convex, it can be shown that
$\|\chi_{\mathcal{D}}(\mathbf{x})-\mathbf{x}_{*}\|\leq\|\mathbf{x}-\mathbf{x}_{*}\|_{2}.$
In many optimization problems, the search space, $\mathcal{D}$, is usually a
$d$-dimensional hyper-rectangle, it follows that $\mathcal{D}$ is convex.

A second step is to enhance an algorithm’s exploration. We achieve this by
injecting noise into its dynamics. In stochastic optimization algorithms, such
as the perturbed gradient descent and the Langevin algorithm, Gaussian noise
is added to the classical gradient descent algorithm. This injection of noise
effectively helps the algorithm to get out of saddle points or local minimums
efficiently [12, 13]. We adopt the same strategy in this article due to its
simplicity. Other form of noise injection also exists. For example, in multi-
arm bandit problem, one can enforce exploration in the classical greedy
algorithm by allocating a small proportion of trials on random arms.

By combining these two steps, the modified algorithm can be formulated as

$\mathbf{x}_{i}(t+1)=\chi_{\mathcal{D}}(\chi_{\mathcal{D}}(\mathbf{x}^{\prime}_{i}(t+1))+\mathbf{w}_{i}(t)),$
(1)

where $\mathbf{x}^{\prime}_{i}(t+1)$ is the position of the $i$th agent after
applying the algorithm’s update rule at $t+1$, and $\mathbf{w}_{i}(t)$ is an
independent draw from a density $p_{t}$ such that there is a constant
$\alpha>0$

$p_{t}(\mathbf{y}-\mathbf{x})>\alpha,\quad\mbox{for
}\mathbf{x},\mathbf{y}\in\mathcal{D}.$ (2)

The noise density $p_{t}$ in general can be adaptive, depending on the
evolution of the algorithm up to time $t$. Here we only need it to be strictly
positive and bounded away from zero. We call this method of modification in
(1) a perturbation-projection (PP) modification or strategy. The proposition
below shows that with PP strategy, the modified algorithm satisfies (C1) and
(C2).

######  Proposition II.2.

Suppose $A$ is a swarm-based algorithm and at each iteration, every agent is
projected onto $\mathcal{D}$. If the perturbation-projection modification (1)
is applied to at least one of the agents, then (C1) and (C2) hold.

###### Proof.

Since $\chi_{\mathcal{D}}$ is applied to all agents, (C1) clearly holds. At
the $(t+1)$th iteration, denote the agent to which the modification is applied
by $i=i(t+1)$. Let
$\mathbf{y}=\chi_{\mathcal{D}}(\mathbf{x}^{\prime}_{i}(t+1))$.

Then, for any $d$-dimensional ball $B$ in $\mathcal{D}$,

$\displaystyle\max_{1\leq k\leq n}\mathbb{P}_{t}(\mathbf{x}_{k}(t+1)\in B)$
$\displaystyle\geq$ $\displaystyle\mathbb{P}_{t}(\mathbf{x}_{i(t+1)}(t+1)\in
B)$ $\displaystyle\geq$
$\displaystyle\mathbb{P}_{t}(\mathbf{y}+\mathbf{w}_{i}(t)\in B)$
$\displaystyle=$
$\displaystyle\int_{B}\mathbb{P}_{t}(\mathbf{y}+\mathbf{w}_{i}(t)\in
d\mathbf{x})$ $\displaystyle=$
$\displaystyle\int_{B}\mathbb{P}_{t}(\mathbf{w}_{i}(t)\in
d(\mathbf{x}-\mathbf{y}))$ $\displaystyle=$
$\displaystyle\int_{B}p_{t}(\mathbf{y}-\mathbf{x})d\mathbf{x}\geq\alpha|B|,$

showing (C2) holds. ∎

##  III Applications

In this section, we illustrate how to implement our PP strategy to three
nature-inspired swarm-based optimization algorithms, namely, particle swarm
optimization (PSO), bat algorithm (BAT) and competitive swarm optimization
(CSO).

###  III-A PSO and mPSO

Kennedy and Eberhart [14] proposed the particle swarm optimization (PSO)
algorithm. The algorithm has stimulated many refinements and inspired many
variants, and these algorithms find wide applications in many fields. A recent
search of “particle swarm optimization” in the Web of Science generated more
than to 32,000 articles and almost 600 review articles. A small sampling of
early examples of how PSO has been studied and modified over the years are
[15, 16, 17, 18, 19, 20, 19, 21, 22, 23]. PSO probably has the most modified
versions among nature-inspired metaheuristic algorithms and they continue to
this day.

PSO algorithm models after the movement of a flock of birds looking for food
collectively. To align with our terminology in this article, we think of birds
in PSO as agents, and food as a global minimizer of an objective function. At
the next iteration, $t+1$, each agent veers towards the best location it has
found (cognitive/memory component) and the best location known to the flock up
to iteration $t$ (social component). To describe the memory effect, we
allocate a personal “memory” agent $\mathbf{x}^{*}_{i}(t)$ to record the best
location agent $i$ has found up to iteration $t$, and a global “memory” agent
$\mathbf{x}^{*}(t)$ to record the best location found by all agents
collectively so far. Recall the basic PSO algorithm runs the following
iterations iteratively:

  1. 1.

Generate two independent random vector $\mathbf{U}_{1}$ and $\mathbf{U}_{2}$
from the uniform distribution on $[0,1]^{d}$ and let

$\displaystyle\mathbf{v}_{i}(t+1)$ $\displaystyle=$ $\displaystyle
w\mathbf{v}_{i}(t)+c_{1}\mathbf{U}_{1}\circ(\mathbf{x}_{i}^{*}(t)-\mathbf{x}_{i}(t))$
$\displaystyle\quad+c_{2}\mathbf{U}_{2}\circ(\mathbf{x}^{*}(t)-\mathbf{x}_{i}(t)),$

where $w,c_{1},c_{2}$ are given constants, and $\circ$ denotes the Hadamard
product.

  2. 2.

Update $\mathbf{x}_{i}(t+1)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)$.

  3. 3.

Update personal best

$\mathbf{x}^{*}_{i}(t+1)=\begin{cases}\mathbf{x}_{i}(t+1),&\mbox{if
}f(\mathbf{x}^{*}_{i}(t))>f(\mathbf{x}_{i}(t+1));\\\
\mathbf{x}^{*}_{i}(t),&\mbox{otherwise}.\end{cases}$

  4. 4.

Repeat steps 1–3 for all agents, then update the global best agent

$\mathbf{x}^{*}(t+1)=\begin{cases}\mathbf{x}_{j}^{*}(t+1),&\mbox{if condition
H holds},\\\ \mathbf{x}^{*}(t),&\mbox{otherwise},\end{cases}$

where $j=\text{argmin}_{1\leq i\leq n}f(\mathbf{x}^{*}_{i}(t+1))$,
$m(t)=\min_{1\leq i\leq n}f(\mathbf{x}^{*}_{i}(t))$ and condition H:
$m(t+1)<m(t)$.

Our modified PSO, denoted by mPSO, is to apply the PP strategy in step 2:

2’) Update

$\mathbf{x}_{i}(t+1)=\chi_{\mathcal{D}}\left(\chi_{\mathcal{D}}\left(\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)\right)+\mathbf{w}_{i}(t)\right),$
(3)

where $\mathbf{w}_{i}$’s are independent multivariate normal distributed with
mean vector $\mathbf{0}$ and covariance matrix $\sigma\mathbf{I}$.

Since its introduction in 1995, there have been numerous attempts to analyze
the basic PSO convergence behavior. In our opinion, Yuan and Yin [24] provide
the first rigorous proof of the weak convergence of PSO to a global optimum,
without overly restrictive assumptions. Moreover, using stochastic
approximation technique, they provide the rate of convergence. In another
direction, Tong et al. [25] introduce two smoothed versions of PSO and prove
that they converge a.s. to a cost function’s global optimum.

The next proposition shows that if we modify the basic PSO algorithm by
replacing the update step 2 by the new noise enhanced update step 2’, the
modified mPSO algorithm converges to a global minimum of the cost function
a.s.

######  Proposition III.1.

Let $f$ be a continuous function defined on a compact subset $\mathcal{D}$ of
$\mathbb{R}^{d}$ and one of its minimizers is in the interior of
$\mathcal{D}$. Then the algorithm mPSO converges to
$\min_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})$ a.s..

###### Proof.

By Theorem II.1, it suffices to verify that mPSO satisfies (C1)–(C3). Note
that $\mathbf{x}_{i}(t+1)\in\mathcal{D}$ by projection. As
$\mathbf{x}_{i}^{*}(t+1)$ takes either value $\mathbf{x}_{i}^{*}(t)$ or
$\mathbf{x}_{i}(t+1)$, therefore $\mathbf{x}_{i}^{*}(t+1)\in\mathcal{D}$ if
$\mathbf{x}_{i}^{*}(t)\in\mathcal{D}$. Similar argument applies to
$\mathbf{x}^{*}(t+1)$, so (C1) holds by induction. Since we add noise to all
agents except the $\mathbf{x}_{i}^{*}(1\leq i\leq n)$ and $\mathbf{x}^{*}$
agents, Proposition II.2 implies (C2) holds. Consider the global memory agent,
by step 4, $f(\mathbf{x}^{*}(t+1))\leq f(\mathbf{x}^{*}(t))$, so (C3) is
verified. This completes the proof. ∎

###  III-B BAT and mBAT

The success of bats, dolphins, shrews and other animals in applying
echolocation technique to hunt and navigate inspired Yang [26] to propose the
Bat Algorithm (BAT) in 2010. BAT meets with rising popularity in applications.
According to a recent search of the Web of Science, over 1,500 articles,
nearly 800 proceedings papers and 48 review papers have been written on this
algorithm and its variants in just ten years. We refer interested readers in
BAT and its wide ranging applications to the review articles in [27, 28] and
the references therein. Recent applications of the BAT algorithm include [29,
30]. As with other metaheuristic algorithms, BAT has also many modified or
hybridized versions for improved performance in various ways; some examples of
the modified versions can be found in [31, 32, 33, 34, 35] and some examples
of BAT algorithm hybridized with another metaheuristic algorithm are [36, 37,
38].

Let $\mathbf{x}_{i}(t)$, $1\leq i\leq n$, denotes the position of the $i$-th
bat at $t$-th iteration, and $\mathbf{x}^{*}(t)$ tracking the best position
found by the cauldron of $n$ bats. Recall BAT runs the following iterations
iteratively:

  1. 1.

Generate $n$ independent random number $U_{i}$ from the uniform distribution
on $[f_{\min},f_{\max}]$,

$\displaystyle\mathbf{v}_{i}(t+1)$ $\displaystyle=$
$\displaystyle\mathbf{v}_{i}(t)+U_{i}(\mathbf{x}_{i}(t)-\mathbf{x}^{*}(t)).$

  2. 2.

Let $r_{0}$ (pulse rate) denote a user-specified threshold. Generate $n$
independent random numbers $r_{i}$ from the uniform distribution on $[0,1]$,
$1\leq i\leq n$.

Update according to (a): $r_{i}<r_{0}$ or (b)$r_{i}\geq r_{0}$:

(a) If $r_{i}<r_{0}$,
$\mathbf{x}_{i}(t+1)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)$;

(b) If $r_{i}\geq r_{0}$,
$\mathbf{x}_{i}(t+1)=\mathbf{x}^{*}(t)+\epsilon_{i}(t)$, where the components
of $\epsilon_{i}(t)$ are independent mean 0 normal distribution with standard
deviation 0.001.

  3. 3.

Generate $n$ independent random numbers $r_{i}$ from the uniform distribution
on $[0,1]$. If $r_{i}$ is less than a threshold $r_{A}$ (loudness), or
$f(\chi_{\mathcal{D}}(\mathbf{x}_{i}(t)))<f(\mathbf{x}_{i}(t+1))$ update
$\mathbf{x}_{i}(t+1)=\chi_{\mathcal{D}}(\mathbf{x}_{i}(t))$.

  4. 4.

Update global best

$\mathbf{x}^{*}(t+1)=\mathbf{x}_{i^{*}}(t+1),$

where $i^{*}=\arg\min_{i}f(\mathbf{x}_{i}(t+1)).$

In the procedure above, $f_{\min},f_{\max}$ are given constants which describe
the minimal and maximal frequencies. The threshold probabilities $r_{0}$ and
$r_{A}$ describe the pulse and emission rates. The random noise
$\epsilon_{i}(t)$ describes the loudness of each bat. Yang [26] also suggests
using time varying $\epsilon_{i}(t)$ and $r_{0}$ for improved performance.
Here we fix them as in the standard values suggested by the MATLAB package
provided by the same author.

Our modified BAT, denoted by mBAT, is to apply the PP strategy by an
additional step after step 3

3’) Update

$\mathbf{x}_{i}(t+1)=\chi_{\mathcal{D}}\left(\chi_{\mathcal{D}}\left(\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)\right)+\mathbf{w}_{i}(t+1)\right).$
(4)

######  Proposition III.2.

Let $f$ be a continuous function defined on $\mathcal{D}$, a compact subset of
$\mathbb{R}^{d}$. Suppose further that a minimizer of $f$ lies in the interior
of $\mathcal{D}$, then the algorithm mBAT converges to
$\min_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})$ a.s..

###### Proof.

By Theorem II.1, it suffices to verify mBAT satisfies (C1)–(C3). Note that
$\mathbf{x}_{i}(t+1)\in\mathcal{D}$ by projection. As $\mathbf{x}^{*}(t+1)$
takes either value $\mathbf{x}^{*}(t)$ or $\mathbf{x}_{i}(t+1)$ for some $i$,
therefore it is in $\mathcal{D}$ if $\mathbf{x}_{i}^{*}(t)\in\mathcal{D}$. So
(C1) holds by induction.

To verify (C2), we denote the value of $\mathbf{x}_{i}(t+1)$ after step 3’)
(4) as $\mathbf{y}_{i}(t)$. Following the proof of Proposition II.2, we can
show that for any d-ball $B\subset\mathcal{D}$,

$\mathbb{P}(\mathbf{y}_{i}(t)\in B)\geq\alpha|B|.$

Then if all $\mathbf{y}\in B$ satisfies $f(\mathbf{y})<f(\mathbf{x}_{i}(t))$,
we find that

$\mathbb{P}(\mathbf{x}_{i}(t+1)\in
B)\geq(1-r_{A})\mathbb{P}(\mathbf{y}_{i}(t)\in B)\geq\alpha(1-r_{A})|B|.$

So (C2) holds.

By step 4, $f(\mathbf{x}^{*}(t+1))\leq f(\mathbf{x}^{*}(t))$, so (C3) is
verified. This completes the proof of Proposition III.2. ∎

###  III-C CSO and mCSO

Inspired by the particle swarm optimization, Cheng and Jin [39] in 2014
proposed the Competitive Swarm Optimizer (CSO) for “large scale optimization
problems and is able to effectively solve problems of dimensionality up to
5000”. It is generally deemed to outperform PSO and it is also continually
improved in various ways such as [40, 41]. Most recently, [42] applied CSO to
optimize a multitasking problem and Zhang et al. [43] modified CSO to find
optimal experiment designs to estimate interesting parameters in a nonlinear
regression model with several interacting factors.

We shall call the particles in CSO agents. In each CSO iteration, agents are
randomly paired and their functional values are compared. In each pair, agent
with smaller functional value is declared winner and its configuration remains
unchanged in the next iteration. The other agent in the pair is declared loser
and his configuration in the next iteration will move closer to the winner’s
configuration and the average configuration of all the agents, see steps 4 and
5 below. Specifically, CSO repeats steps 1 to 5 until termination criterion is
met.

  1. 1.

Randomly pair the agents. Repeat the steps 2 to 5 below for each pair.

  2. 2.

Consider a pair of agents indexed as $i$ and $j$. Define $w=i,\ell=j$ if
$f(\mathbf{x}_{i}(t))<f(\mathbf{x}_{j}(t))$. Otherwise $w=j,\ell=i$.

  3. 3.

Let $\mathbf{x}_{w}(t+1)=\mathbf{x}_{w}(t)$,
$\mathbf{v}_{w}(t+1)=\mathbf{v}_{w}(t)$.

  4. 4.

Generate three independent random vectors $U_{1},U_{2},U_{3}$ from the uniform
distribution on $[0,1]^{d}$ and let

$\displaystyle\mathbf{v}_{\ell}(t+1)$ $\displaystyle=$ $\displaystyle
U_{1}\circ\mathbf{v}_{\ell}(t)+U_{2}\circ(\mathbf{x}_{w}(t)-\mathbf{x}_{\ell}(t))$
$\displaystyle\quad+\phi U_{3}\circ(\bar{\mathbf{x}}(t)-\mathbf{x}_{l}(t)),$

where $\circ$ denotes element-wise product, $\phi$ a given constant, and
$\bar{\mathbf{x}}(t)=\sum_{i=1}^{n}\mathbf{x}_{i}(t)/n$.

  5. 5.

Update $\mathbf{x}_{\ell}(t+1)=\mathbf{x}_{\ell}(t)+\mathbf{v}_{\ell}(t+1)$.

Our modified CSO, denoted by mCSO, is to apply the PP strategy in step 5 only
to the loser agent:

5’) Update

$\mathbf{x}_{\ell}(t+1)=\chi_{\mathcal{D}}\left(\chi_{\mathcal{D}}\left(\mathbf{x}_{\ell}(t)+\mathbf{v}_{\ell}(t+1)\right)+\mathbf{w}_{\ell}(t)\right).$
(5)

It is not necessary to apply projection on $\mathbf{x}_{w}(t+1)$ because
$\mathbf{x}_{w}(t+1)=\mathbf{x}_{w}(t)\in\mathcal{D}$.

######  Proposition III.3.

Let $f$ be a continuous function over $\mathcal{D}$, a compact subset of
$\mathbb{R}^{d}$. Suppose further that a minimizer of $f$ lies in the interior
of $\mathcal{D}$, then the algorithm mCSO converges to
$\min_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})$ a.s..

###### Proof.

By Theorem II.1 and Proposition II.2, our claim will be proved if we can show
that mCSO satisfies (C1)–(C3). Since $\mathbf{x}_{\ell}(t+1)\in\mathcal{D}$ by
projection, and $\mathbf{x}_{w}(t+1)=\mathbf{x}_{w}(t)$, (C1) holds. As we add
noise to the loser agents, (C2) holds by Proposition II.2. Let
$i_{*}=\arg\min_{1\leq i\leq n}f(\mathbf{x}_{i}(t))$. In the next iteration,
$\mathbf{x}_{i_{*}}(t)$ will be a winner no matter whom he is paired with. So
$\mathbf{x}_{i_{*}}(t+1)=\mathbf{x}_{i_{*}}(t)$; and
$\min_{i}f(\mathbf{x}_{i}(t+1))\leq
f(\mathbf{x}_{i_{*}}(t+1))=\min_{i}f(\mathbf{x}_{i}(t)).$ This proves that
(C3) is satisfied and hence the proof of the proposition. ∎

##  IV Heterogeneous PP

In applications, the strength of stochastic perturbation in the PP strategy
often needs proper tuning. In general, a stronger perturbation enhances
exploration at the expense of exploitation. During the final stage of
optimization, the agents will have difficulties in narrowing down the optimal
solution with strong perturbations. There can be different ways to solve this
issue such as choosing a diminishing perturbation as in the original BAT
algorithm; using adaptive adjustment mechanism; or using reinforcement
learning tools to find the optimal parameter. However, these approaches lead
to further tuning parameters or additional computational complexity. Instead,
we propose a simpler strategy, which is nature-inspired, that involves role
heterogeneity of members in the swarm.

For swarms in nature, it is common to observe they adoption division of
labors. Agents specializing in different tasks collaborate to improve
collective effectiveness. In the same vein, we set some agents in the swarm to
specialize in exploration while others in exploitation. In the context of PP
strategy, we only perturb the exploration agents so they help exploring the
search space while we do not perturb exploitation agents so as not to hamper
their convergence capability. We call this heterogeneous perturbation-
projection, abbreviated to HPP, strategy. In contrast with the other ways to
find optimal perturbation strength, HPP strategy is very easy to implement.
For example, if we want to have half exploration agents and half exploitation
agents, we will apply PP strategy only to the first half of the agents while
the other agents we only apply the projection modification. This idea can be
applied directly to mPSO and mBAT. As for mCSO, because pairs are randomly
formed, we will apply PP strategy only to the loser agents in the first half
of pairs.

Our HPP strategy can be considered as a special case of cooperative or
collaborative learning, which suggests that mixing agents of different
settings can lead to improved overall performance. Existing works on
cooperative swarm-based algorithms can be found in [20, 44, 45]. In
comparison, our HPP strategy focuses more on the cooperation between
exploration agents and exploitation agents. Recently, [46] also applies this
idea to local optimization algorithms such as gradient descent and Langevin
algorithm.

Note that Proposition II.2 only requires the perturbation strategy applied to
at least one agent, so this proposition is still applicable to algorithms with
HPP, which perturbs half of their agents (or loser agents in CSO).
Consequently, Theorem II.1 guarantees H-PP modified algorithms converge a.s to
a global optimal solution. We denote the modified version of an algorithm $A$
using H-PP strategy by $hmA$.

##  V Numerical experiments

In this section, we conduct extensive numerical experiments to investigate
whether our modified versions of $A$ ($A=$ PSO, BAT or CSO), $mA$ and $hmA$,
outperform or at least perform on par with $A$ given a computing budget for a
large class of very different types of cost functions. We use two measures to
evaluate the performance of the modified algorithms relative to the original
algorithm. They measure how likely and by how much the modified algorithm
obtains better results than the original algorithm.

###  V-A Test functions and details of experiments

We use a total of 30 test functions of varying dimensions listed in [47, Table
2] for our experiments. The list of test functions was compiled from various
sources by the authors in [47]. They include commonly used test functions such
as Ackley, Griewank, Powell, Rastrigin, Sphere etc; and are of different
kinds, unimodal, multi-modal, separable and inseparable. Majority of the
functions in their table can be defined for arbitrary dimension $d$. For such
functions, we consider $d=5,10,20$ and $40$ to investigate the effect of
dimension on the optimization algorithms. The other test functions (except one
which can only be defined for dimension 4, which we exclude in our
experiments) can only be defined for $d=2$. In total, we have a collection of
70 test functions, see Appendix A, denoted by $\mathcal{C}$. We shall
partition $\mathcal{C}$ into different collections according to $d=2,5,10,20$
and $40$ in reporting our results.

For each $f\in\mathcal{C}$, we conduct 100 runs of algorithms $A$, $mA$ and
$hmA$ for 10,000 iterations where $A$ = PSO, BAT, or CSO. To gauge the
progress of each run of an algorithm, we output the algorithm best functional
value found at $t$-th iteration, where $t=50,100,200,400,1000,3000$ and
$10000$.

Throughout the experiments, the number of agents/particles used is $n=32$. For
$mA$ or $hmA$ ($A=$ PSO, BAT or CSO), the stochastic perturbation used is
independent normal distributions with mean 0 and standard deviation $0.005$.

We follow common parameter settings for PSO, BAT and CSO:

  * •

PSO: $w=0.729$ (inertia weight), $c_{1}=c_{2}=1.5$ (acceleration constants);

  * •

BAT: $Q_{\min}=0$ (minimum frequency), $Q_{\max}=100$ (maximum frequency),
$r_{0}=0.5$ (pulse rate), $r_{A}=0.5$ (loudness);

  * •

CSO: $\phi=0$.

###  V-B Comparison and results

This subsection compares the performance of an algorithm $A$ versus one of its
modification $B$ using several test functions taken from the set
$\mathcal{F}$. The functions in the set are commonly used to compare
performance of an algorithm in the engineering literature. To study how the
performance on a problem depends on the dimension of the problem, each
$\mathcal{F}$ here is the subset of $\mathcal{C}$ that contains functions of
the same dimension $d=2,5,10,20$ or $40$. We are interested in the following
two questions: How likely does $B$ outperform $A$? How much does $B$ exceed
the best possible value from using $A$ and $B$? Here $B=mA$ or $hmA$.

(a) How likely will $B$ outperform $A$?

Let $A_{f}(t)$ and $B_{f}(t)$ denote the (theoretical) best functional values
from algorithms $A$ and $B$ respectively for the test function $f$ up to the
$t$-th iteration. Similarly, let $A_{r,f}(t)$ and $B_{r,f}(t)$ denote the best
functional values output by algorithms $A$ and $B$ respectively up to the
$t$-th iteration at the $r$-th run. The probability that algorithm $B$
outperforms algorithm $A$ for $f$ at the $t$-th iteration is
$\mathbb{P}(B_{f}(t)<A_{f}(t))$, and it can be estimated by
$\sum_{r=1}^{100}I(B_{r,f}(t)<A_{r,f}(t))/100$. The average of these estimates
over $\mathcal{F}$, denoted by $P_{B\succ A}(t)$, can be thought of as an
average winning proportion of $B$ over $A$ (or simply, it’s winning
proportion). Specifically,

$P_{B\succ
A}(t):=\frac{1}{100|\mathcal{F}|}\sum_{f\in\mathcal{F}}\sum_{r=1}^{100}I(B_{r,f}(t)<A_{r,f}(t)).$

Similarly, interchanging $A$ and $B$ above, we define the winning proportion
of $A$ over $B$, $P_{A\succ B}(t)$. Since $P_{A\succ B}(t)+P_{B\succ A}(t)=1$,
it suffices to track $P_{B\succ A}(t)$. Loosely speaking, if $P_{B\succ
A}(t)>1/2$, $B$ is more likely to outperform $A$ for $f\in\mathcal{F}$.
Indeed, the larger the $P_{B\succ A}(t)$, the more likely $B$ outperforms $A$.

Figure 1: Plots of winning proportion of $mA$ against $A$ (top row), and that
of $hmA$ against $A$ (bottom row) where $A=$ PSO, BAT and CSO categorized
according to the dimension of the test functions.

The numerical results are plotted in Figure 1. Based on these plots, we
observe have the followings. When $t$ is fixed, the winning proportion of
PP/HPP modified algorithm outperforming its counterpart increases as the
dimension of the test function increases. Both mBAT and hmBAT are comparable
in their performance against BAT, and the winning proportion increases
steadily from about 0.5 to 0.7 for dimension at least 5. Interestingly, the
HPP strategy has more significant enhancement effect on PSO and CSO. For test
functions with dimension 40, winning proportion of hmCSO is nearly 1 after
3000 iterations.

(b) How much does $B$ (or $A$) outperform $A$ (or $B$)?

For simplicity, we suppress the dependence of a notation on other
parameters/notations when the context is clear. To this end, let $m_{*}$ and
$m^{*}$ denote respectively, the minimum (i.e., the best) and maximum (i.e.,
the worst) of the outputs from 100 runs of $A$ and 100 runs of $B$ for $f$ at
the $t$-th iteration. We view $m_{*}$ as a proxy of the best possible output
by $A$ and $B$; and let $R=m^{*}-m_{*}$ denote the range of the 200 outputs of
$A$ and $B$. Define the relative error of $A$ with respect to $A$ and $B$,
$RE_{A,A\cdot B,f}(t)$, to be the average of $(A_{r,f}(t)-m_{*})/R$ over 100
runs; and the relative error of $B$ relative to $A$ and $B$, $RE_{B,A\cdot
B,f}(t)$ is similarly defined. We consider the relative instead of absolute
error in order to cancel out the effect due to a scale change of $f$ (i.e., if
we consider $cf$ instead of $f$). For brevity, we simply call this quantity
the relative error. Note that the relative errors lie between 0 and 1. Small
relative error of $A$ (i.e., $RE_{A,A\cdot B,f}(t)$) indicates the results
from algorithm $A$, as a whole, are closer to the best possible value $m_{*}$.
Further, if the relative error of $A$ is less than that of $B$, then generally
results from algorithm $A$ are closer to $m_{*}$ than those by $B$’s. A
similar interpretation can be extended over a class of functions
$f\in\mathcal{C}$ if we define $RE(A,A\cdot
B)(t):=\frac{1}{|\mathcal{F}|}\sum_{f\in\mathcal{F}}RE_{A,A\cdot B,f}(t)$ and
call this the overshoot of $A$ relative to $A$ and $B$.

Figure 2: Plots of relative error of $A$ relative to $A$ and $mA$ (dashed blue
curve), and relative error of $mA$ relative to $A$ and $mA$ (solid blue
curves). Plots of relative error of $A$ relative to $A$ and $hmA$ (dashed red
curves), and relative error of $hmA$ relative to $A$ and $hmA$ (solid red
curves).

Figure 2 plots the relative errors of $A$ (dashed blue curve) and $mA$ (solid
blue curve) with respect to one another; and the relative errors of $A$
(dashed red curve) and $hmA$ (solid red curve) with respect to one another. We
only display the plots for $d=10$ and $d=40$. We highlight some observations
from these plots in Figure 2. (i) hmBAT and mBAT outperform BAT by a
significant margin. (ii) At $d=40$, both mPSO and hmPSO are nearly 0 from 1000
iterations onwards. This implies mPSO and hmPSO significantly outperform the
basic PSO from 1000 iterations. At $d=10$, only hmPSO (resp., hmCSO) improves
the basic PSO (resp., CSO) from 400 iterations.

Combining the observations in (a) and (b), Figures 1 and 2 convincingly
demonstrate that the HPP strategy significantly improves BAT, CSO and PSO
particularly when the dimension of a test function is large. Improvement is
measured in two complementary aspects, “how likely” and “how much” one
algorithm outperforms the other.

###  V-C Tuning parameter analysis for PP and HPP strategies

As demonstrated in Section V-B, additional stochastic component in a modified
swarm-based algorithm enhances the exploration ability of the algorithm under
consideration. The stochastic component chosen is Gaussian with mean 0 and
standard deviation, $\sigma=0.005$. Expectedly the noise level, in terms of
$\sigma$ or a heavier-tailed distribution, affects the efficiency of the
modified algorithm since strong level of noise interferes the algorithm’s
exploitation ability. We conduct further numerical experiments with the same
set-up as above to investigate the effects of $\sigma$ in the Gaussian noise
and $t$-distributions, which have heavier tails than the Gaussian, on the
performance of $mA$ and $hmA$.

(i) On the standard deviation of the Gaussian noise We tried
$\sigma=0.005,0.01,0.02$ and $0.05$. The plots of $P_{B\prec A}(t)$ (analogous
to Figure 1) and $RE(A,A\cdot B)(t)$ (analogous to Figure 2) are very similar
for $hmA$ ($A=$ PSO, BAT and CSO) for all the $\sigma$s; and $\sigma=0.005$
and $0.01$. Hence, we recommend using $\sigma$ in $[0.005,0.01]$. Interested
readers can refer to Figures S4 and S5 in Appendix B.

(ii) On the choice of distribution We examine the effect on the performance of
the modified algorithms if we change the Gaussian distribution with a heavier-
tailed distribution. We choose distributions $0.01\sqrt{(m-2)/m}\ t_{m}$ for
$m=5,10,30$ and $60$ where $t_{m}$ denotes a $t$-distribution with $m$ degrees
of freedom. The scaling factor $0.01\sqrt{(m-2)/m}$ is to ensure the scaled
$t$-distribution has the same standard deviation as that of the Gaussian noise
in (i) above. Note also that $t_{m}$ converges to a standard normal. For $m$
large, we do not expect essential difference between the $t$-distribution and
the Gaussian noise.

The analogous plots of $P_{B\prec A}(t)$ and $RE(A,A\cdot B)(t)$ using
$t$-distributions for the choices of degrees of freedom (df) look almost the
same as the corresponding plots using Gaussian noise. Interested readers can
view these plots as Figure S6-S7 in Appendix B.

##  VI Conclusion

In this article, we delineate sufficient conditions for a swarm-based
algorithm to converge almost surely to an optimum of an objective function.
This assumes that there is no constraint on the computation budget and the
algorithm is allowed to run indefinitely. The class of objective functions for
which Theorem II.1 holds is very large. We proposed two strategies,
perturbation-project (PP) and heterogeneous perturbation-projection (HPP), to
modify a given swarm-based algorithm so that the sufficient conditions
(C1)–(C3) hold, and hence Theorem II.1 guarantees the modified algorithm
converges almost surely to the global optimum. The proposed strategies are
simple to implement and not algorithm-specific. We demonstrated how these
strategies are applied to PSO, Bat Algorithm and CSO. Extensive numerical
experiments were conducted and the results show that the modified algorithm
either out-performs or performs on par with the original algorithm over a
finite computational budget, and especially so for HPP modified algorithm. Our
simulation results suggest that applying the HPP strategy with mean 0 and
standard deviation 0.005 Gaussian distribution for the additional stochastic
component is likely to be effective in improving a swarm-based algorithm.

To conclude this article, we report two interesting observations from our
experiments. The first observation is that when $A$ outperforms $mA$ or $hmA$,
$A$ frequently produces a marginally better result; however, when $mA$ or
$hmA$ outperforms $A$, the modified version quite often outperforms $A$ by
relatively larger margin. We offer a heuristic argument. When both $A$ and
$mA$ are exploring in the same neighbourhood of a local minimum, the
stochastic element hinders the exploitation by a small amount and hence $A$ is
more likely to produce only marginally better solution. However, when the
occasion for $mA$ to make a big jump leaving the local minimum’s neighbourhood
arises, $mA$ has the potential to explore a better region and hence produces
noticeably better solution than the original algorithm does. The same
heuristic reasoning carries over to the comparison of $A$ and $hmA$. Moreover,
$hmA$ version due to the built in diversity of agents, better balances the
conflicting demands of exploitation and exploration, and expectedly amplifies
this observation.

The second observation comes from our application of PP and HPP strategies to
Differential Evolution (DE), an evolutionary algorithm. The same numerical
experiments were conducted to compare DE with mDE, and DE with hmDE. We found
that DE outperforms mDE and hmDE; the plots can be found in Appendix C as
Figures S8-S9. A reason may be that DE already performs very well in the
optimization problem over majority of the test functions in the experiment and
so any additional stochastic component will either not help or even worsen the
search solution. Moreover, the random perturbation we apply here mimic the
random local exploration of an ant or other animal. Genetic mutation in
general cannot describe such dynamics and require some other more structured
mechanism.

## Acknowledgment

Choi’s research was supported by the Singapore MOE Academic Research Funds
R-155-000-222-114. Lai’s research was supported by the National Science
Foundation under DMS-1811818. Tong’s research was supported by MOE Academic
Research Funds R-146-000-292-114. Wong’s research was supported by the
National Institutes of Health under R01GM107639.

## References

  * Whitacre [2011a] J. M. Whitacre, “Recent trends indicate rapid growth of nature-inspired optimization in academia and industry,” _Computing_ , vol. 93, pp. 121–133, 2011a.
  * Whitacre [2011b] ——, “Survival of the flexible: explaining the recent popularity of nature-inspired optimization within a rapidly evolving world,” _Computing_ , vol. 93, pp. 135–146, 2011b.
  * Huang et al. [2019] C. Huang, Y. Li, and X. Yao, “A survey of automatic parameter tuning methods for metaheuristics,” _IEEE Transactions on Evolutionary Computation_ , vol. 24, no. 2, pp. 201–216, 2019.
  * Choi et al. [2021] K. P. Choi, T. L. Lai, X. Tong, and W. K. Wong, “A statistical approach to adaptive parameter tuning in nature-inspired optimization and optimal sequential design of dose-finding trials,” _Statistical Sinica_ , p. to appear, 2021.
  * Parsopoulos et al. [2001] K. E. Parsopoulos, V. Plagianakos, G. Magoulas, and M. Vrahatis, “Objective function” stretching” to alleviate convergence to local minima,” _Nonlinear Analysis-Theory Methods and Applications_ , vol. 47, no. 5, pp. 3419–3424, 2001.
  * Stacey et al. [2003] A. Stacey, M. Jancic, and I. Grundy, “Particle swarm optimization with mutation,” in _The 2003 Congress on Evolutionary Computation, 2003. CEC’03._ , vol. 2. IEEE, 2003, pp. 1425–1430.
  * Krohling [2005] R. A. Krohling, “Gaussian particle swarm with jumps,” in _2005 IEEE Congress on Evolutionary Computation_ , vol. 2. IEEE, 2005, pp. 1226–1231.
  * Elshamy et al. [2007] W. Elshamy, H. M. Emara, and A. Bahgat, “Clubs-based particle swarm optimization,” in _2007 IEEE Swarm Intelligence Symposium_. IEEE, 2007, pp. 289–296.
  * James et al. [2013] G. James, D. Witten, T. Hastie, and R. Tibshirani, _An introduction to statistical learning_. Springer, 2013, vol. 112.
  * Shalev-Shwartz and Ben-David [2014] S. Shalev-Shwartz and S. Ben-David, _Understanding machine learning: From theory to algorithms_. Cambridge University Press, 2014.
  * Törn and Žilinskas [1989] A. Törn and A. Žilinskas, _Global optimization_. Springer, 1989, vol. 350.
  * Ge et al. [2015] R. Ge, F. Huang, C. Jin, and Y. Yuan, “Escaping from saddle points—online stochastic gradient for tensor decomposition,” in _Proceedings of the Conference on Learning Theory_ , 2015.
  * Chen et al. [2020] X. Chen, S. S. Du, and X. T. Tong, “On stationary-point hitting time and ergodicity of stochastic gradient langevin dynamics.” _Journal of Machine Learning Research_ , vol. 21, no. 68, pp. 1–41, 2020.
  * Kennedy and Eberhart [1995] J. Kennedy and R. C. Eberhart, “Particle swarm optimization,” in _Proceedings of the IEEE International Conference on Neural Networks IV_. IEEE, 1995, pp. 1942–1948.
  * Shi and Eberhart [1998] Y. Shi and R. Eberhart, “A modified particle swarm optimizer,” _1998 IEEE International Conference on Evolutionary Computation Proceedings. IEEE World Congress on Computational Intelligence_ , p. 10.1109/ICEC.1998.699146, 1998\.
  * Bratton and Kennedy [2007] D. Bratton and J. Kennedy, “Defining a standard for particle swarm optimization,” in _2007 IEEE swarm intelligence symposium_. IEEE, 2007, pp. 120–127.
  * Poli [2008] R. Poli, “Analysis of the publications on the applications of particle swarm optimization,” _J. Artificial Evolution and Applications_ , vol. 2008, 2008\.
  * Eberhart and Shi [2000] R. C. Eberhart and Y. Shi, “Comparing inertia weights and constriction factors in particle swarm optimization,” in _Proceedings of the IEEE Congress on Evolutionary Computation_. IEEE, 2000, pp. 84–88.
  * Clerc and Kennedy [2002] M. Clerc and J. Kennedy, “The particle swarm–explosion, stability, and convergence in a multidimensional complex space,” _IEEE: Transactions on Evolutionary Computation_ , vol. 6, pp. 58–72, 2002.
  * Van den Bergh and Engelbrecht [2000] F. Van den Bergh and A. P. Engelbrecht, “Cooperative learning in neural networks using particle swarm optimizers,” _South African Computer Journal_ , vol. 2000, no. 26, pp. 84–90, 2000.
  * Trelea [2003] I. C. Trelea, “The particle swarm optimization algorithm: convergence analysis and parameter selection,” _Information Processing Letters_ , vol. 85, no. 6, pp. 317–325, 2003.
  * Chen and Li [2007] X. Chen and Y. Li, “A modified PSO structure resulting in high exploration ability with convergence guaranteed,” _IEEE Trans. Syst. Man Cybern. B Cybern_ , vol. 37, no. 5, pp. 1271–1289, 2007.
  * Pedersen and Chipperfield [2010] M. Pedersen and A. J. Chipperfield, “Simplifying particle swarm optimization,” _Applied Soft Computing_ , vol. 10, no. 2, pp. 618–628, 2010\.
  * Yuan and Yin [2015] Q. Yuan and G. Yin, “Analyzing convergence and rates of convergence of particle swarm optimization algorithms using stochastic approximation methods,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 7, pp. 1760–1773, 2015.
  * Tong et al. [2021] X. Tong, K. P. Choi, T. L. Lai, and W. K. Wong, “Stability bounds and almost sure convergence of improved particle swarm optimization methods,” _Research in Mathematical Sciences_ , p. to appear, 2021.
  * Yang [2010] X.-S. Yang, “A new metaheuristic bat-inspired algorithm,” in _Nature inspired cooperative strategies for optimization (NICSO 2010)_. Springer, 2010, pp. 65–74.
  * Akhtar et al. [2012] S. Akhtar, A. Ahmad, and E. M. Abdel-Rahman, “A metaheuristic bat-inspired algorithm for full body human pose estimation,” in _2012 Ninth Conference on Computer and Robot Vision_. IEEE, 2012, pp. 369–375.
  * Yang and He [2013] X.-S. Yang and X. He, “Bat algorithm: literature review and applications,” _International Journal of Bio-inspired computation_ , vol. 5, no. 3, pp. 141–149, 2013.
  * Cai et al. [2014] X. Cai, W. Li, L. Wang, Q. Kang, Q. Wu, and X. Huang, “Bat algorithm with gaussian walk for directing orbits of chaotic systems,” _International Journal of Computing Science and Mathematics_ , vol. 5, p. Iss. 2. 10.1504/ijcsm.2014.064070, 2014.
  * Xue et al. [2015] F. Xue, Y. Cai, Y. Cao, Z. Cui, and F. Li, “Optimal parameter settings for bat algorithm,” _International Journal of Bio-Inspired Computation_ , p. 7. 125\. 10.1504/IJBIC.2015.069304, 2015.
  * Osaba et al. [2018] E. Osaba, X. Yang, I. J. Fister, J. Del Ser, P. Lopez-Garcia, and A. Vazquez-Pardavila, “A discrete and improved bat algorithm for solving a medical goods distribution problem with pharmacological waste collection,” _Swarm and Evolutionary Computation_ , p. BASE DATA (2018). doi: 10.1016/j.swevo.2018.04.001, 2018.
  * Binu and Selvi [2015] D. Binu and M. Selvi, “Bfc: Bat algorithm based fuzzy classifier for medical data classification,” _Journal of Medical Imaging and Health Informatics_ , vol. 5, pp. 599–606, 2015.
  * Wang et al. [2016] G. Wang, H. Chu, and S. Mirjalili, “Three-dimensional path planning for ucav using an improved bat algorithm,” _Aerospace Science and Technology_ , vol. 49, p. 10.1016/j.ast.2015.11.040, 2016.
  * Khooban and Niknam [2015] M. Khooban and T. Niknam, “A new intelligent online fuzzy tuning approach for multi-area load frequency control: Self adaptive modified bat algorithm,” _International Journal of Electrical Power & Energy Systems_, vol. 71, pp. 254–261, 2015.
  * Lu et al. [2020] S. Lu, K. Xia, and S. Wang, “Diagnosis of cerebral microbleed via vgg and extreme learning machine trained by gaussian map bat algorithm,” _Journal of Ambient Intelligence and Humanized Computing_ , pp. https://doi.org/10.1007/s12 652–020–01 789–3, 2020.
  * He et al. [2014] X. He, W. Ding, and X. Yang, “Bat algorithm based on simulated annealing and gaussian perturbations,” _Neural Comput & Applic_, vol. 25, p. 459–468, 2014.
  * Shrichandran et al. [2017] G. Shrichandran, S. Sathiyamoorthy, P. Malarchelvi, and S. Kezia, “A hybrid glow-worm swarm optimization with bat algorithm based retinal blood vessel segmentation,” _Journal of Computational and Theoretical Nanoscience_ , vol. 14, pp. 2601–2611, 2017.
  * Kishore et al. [2015] P. Kishore, S. Kishore, E. Kumar, K. Kumar, and P. Aparna, “Medical image watermarking with dwt-bat algorithm,” _2015 International Conference on Signal Processing and Communication Engineering System_ , pp. 270–275, 2015.
  * Cheng and Jin [2014] 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.
  * Mohapatra et al. [2017] P. Mohapatra, K. Nath, and S. Roy, “A modified competitive swarm optimizer for large scale optimization problems,” _Applied Soft Computing_ , vol. 59, pp. 340–362, 2017.
  * Wang et al. [2019] B. Wang, Y. Li, F. Yang, and X. Xia, “A competitive swarm optimizer-based technoeconomic optimization with appliance scheduling in domestic pv-battery hybrid systems,” _Complexity_ , p. Article ID 4824837 — https://doi.org/10.1155/2019/4824837, 2019.
  * Li et al. [2020] W. Li, J. Yuan, H. Luo, Z. Lei, and Q. Xu, “Enhanced competitive swarm optimizer for multi-task optimization,” _UCCOR 2929L Proceedings of the 2020 9th International Conference on Computing and Pattern Recognition_ , p. 455–459 ; https://doi.org/10.1145/3436369.3437417, 2020.
  * Zhang et al. [2020] Z. Zhang, W. K. Wong, and K. C. Tan, “Competitive swarm optimizer with mutated agents for finding optimal designs for nonlinear regression models with multiple interacting factors,” _Memetic Computing_ , vol. 12, no. 3, pp. 219–233, 2020.
  * Van den Bergh and Engelbrecht [2004] F. Van den Bergh and A. P. Engelbrecht, “A cooperative approach to particle swarm optimization,” _IEEE Transactions on Evolutionary Computation_ , vol. 8, no. 3, pp. 225–239, 2004.
  * Zhang et al. [2019] X. Zhang, K.-J. Du, Z.-H. Zhan, S. Kwong, T.-L. Gu, and J. Zhang, “Cooperative coevolutionary bare-bones particle swarm optimization with function independent decomposition for large-scale supply chain network design with uncertainties,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 10, pp. 4454–4468, 2019.
  * Dong and Tong [2020] J. Dong and X. T. Tong, “Replica exchange for non-convex optimization,” _arXiv preprint arXiv:2001.08356_ , 2020.
  * Ab Wahab et al. [2015] M. N. Ab Wahab, S. Nefti-Meziani, and A. Atyabi, “A comprehensive review of swarm optimization algorithms,” _PloS One_ , vol. 10, no. 5, p. e0122827, 2015.

##  Appendix A Test functions used in the numerical simulations

TABLE S1: Table of test functions containing their function label, name, known minimum value, unimodality and separability. Functions F4, 5, 6, 7, 9, 10, 11, 13, 14, 17, 18, 20, 21 have dimension 2, and F15 is dimension 5. All other functions were tested at dimensions 5, 10, 20 and 40. Label | Name | min f(x) | Unimodal | Separable
---|---|---|---|---
F1 | Ackley | 0 | No | No
$-20\exp(-0.2\sqrt{\frac{1}{n}\sum_{i=1}^{n}x_{i}^{2}})-\exp(\frac{1}{n}\sum_{i=1}^{n}\cos(2\pi
x_{i}))+20+e$
F2 | Bohachevsky2 | 0 | No | No
$\sum_{i=1}^{d-1}\left[x_{i}^{2}+2x_{i+1}^{2}-0.3\cos(3\pi x_{i})\cos(4\pi
x_{i+1})+0.3\right]$
F3 | Bohachevsky3 | 0 | No | No
$\sum_{i=1}^{d-1}\left[x_{i+1}^{2}+2x_{i+2}^{2}-0.3\cos(3\pi x_{i}+4\pi
x_{i+1})+0.3\right]$
F4 | Bukin6 | 0 | No | No
$100\sqrt{|x_{2}-0.01x_{1}^{2}|}+0.01|x_{1}+10|$
F5 | DropWave | -1 | No | No
$-\frac{1+\cos(12\sqrt{x_{1}^{2}+x_{2}^{2}})}{0.5(x_{1}^{2}+x_{2}^{2})+2}$)
F6 | Eggholder | -959.6407 | No | No
$-(x_{2}+47)\sin(\sqrt{|x_{2}+\frac{x_{1}}{2}+47|})-x_{1}\sin(\sqrt{|x_{1}-(x_{2}+47)|})$)
F7 | GoldSteinPrice | 3 | No | No
$[1+(x_{1}+x_{2}+1)^{2}(19-14x_{1}+3x_{1}^{2}-14x_{2}+6x_{1}x_{2}+3x_{2}^{2})]\times$
$[30+(2x_{1}-3x_{2})^{2}(18-32x_{1}+12x_{1}^{2}+4x_{2}-36x_{1}x_{2}+27x_{2}^{2})]$
F8 | Griewank | 0 | No | No
$1+\sum_{i=1}^{d}\frac{x_{i}^{2}}{4000}-\prod_{i=1}^{d}\cos(\frac{x_{i}}{\sqrt{i}})$
F9 | McCormick | -1.9133 | No | No
$\sin(x_{1}+x_{2})+(x_{1}-x_{2})^{2}-1.5x_{1}+2.5x_{2}+1$
F10 | Schaffer2 | 0 | No | No
$0.5+\frac{\sin^{2}(x_{1}^{2}-x_{2}^{2})-0.5}{(1+0.001(x_{1}^{2}+x_{2}^{2}))^{2}}$
F11 | Schaffer4 | 0.292579 | No | No
$0.5+\frac{\cos^{2}(\sin(\lvert
x_{1}^{2}-x_{2}^{2}\rvert))-0.5}{(1+0.001(x_{1}^{2}+x_{2}^{2}))^{2}}$)
F12 | Bohachevsky1 | 0 | No | Yes
$\sum_{i=1}^{d-1}x_{i}^{2}+2x_{i+1}^{2}-0.3\cos(3\pi x_{i})-0.4\cos(4\pi
x_{i+1})+0.7$
F13 | Booth | 0 | No | Yes
$(x_{1}+2x_{2}-7)^{2}+(2x_{1}+x_{2}-5)^{2}$
F14 | Branin | 0.397887 | No | Yes
$(x_{2}-\frac{5.1}{4\pi^{2}}x_{1}^{2}+\frac{5}{\pi}x_{1}^{2}-6)^{2}+10(1-\frac{1}{8\pi}\cos
x_{1}+10$)
F15 | Michalewicz5 | -4.687658 | No | Yes
$-\sum_{i=1}^{5}\sin(x_{i})\sin^{20}(\frac{ix_{i}^{2}}{\pi})$
F16 | Rastrigin | 0 | No | Yes
$10d+\sum_{i=1}^{d}[x_{i}^{2}-10\cos(2\pi x_{i})]$
F17 | Shubert | -186.73 | No | Yes
$\prod_{j=1}^{2}(\sum_{i=1}^{5}i\cos((i+1)x_{j}+i))$
F18 | Beale | 0 | Yes | No
$(1.5-x_{1}+x_{1}x_{2})^{2}+(2.25-x_{1}+x_{1}x_{2}^{2})^{2}+(2.625-x_{1}+x_{1}x_{2}^{3})^{2}$
F19 | DixonPrice | 0 | Yes | No
$(x_{1}-1)^{2}+\sum_{i=2}^{d}i(2x_{i}^{2}-x_{i-1})^{2}$)
F20 | Easom | -1 | Yes | No
$-\cos(x_{1})\cos(x_{2})\exp(-(x_{1}-\pi)^{2}-(x_{2}-\pi)^{2})$)
F21 | Matyas | 0 | Yes | No
$0.26(x_{1}^{2}+x_{2}^{2})-0.48x_{1}x_{2}$
F22 | Powell | 0 | Yes | No
$\sum_{i=1}^{d/4}[(x_{4i-3}+10x_{4i-2})^{2}+5(x_{4i-1}-x_{4i})^{2}+(x_{4i-2}-2x_{4i-1})^{4}+10(x_{4i-3}-x_{4i})^{4}]$
F23 | Rosenbrock | 0 | Yes | No
$\sum_{i=1}^{d-1}[100(x_{i+1}-x_{i}^{2})^{2}+(x_{i}-1)^{2}]$
F24 | Schwefel | -418.9829d | Yes | No
$418.9829d-\sum_{i=1}^{d}x_{i}\sin(\sqrt{|x_{i}|})$
F25 | Trid6 | -d(d+4)(d-1)/6 | Yes | No
$\sum_{i=1}^{d}(x_{i}-1)^{2}-\sum_{i=1}^{d}x_{i}x_{i-1}$
F26 | Zakharov | 0 | Yes | No
$\sum_{i=1}^{d}x_{i}^{2}+(\sum_{i=1}^{d}0.5ix_{i})^{2}+(\sum_{i=1}^{d}0.5ix_{i})^{4}$
F27 | Sphere | 0 | Yes | Yes
$\sum_{i=1}^{d}x_{i}^{2}$
F28 | Sumsquare | 0 | Yes | Yes
$\sum_{i=1}^{d}ix_{i}^{2}$

##  Appendix B Additional numerical results

Here we report the simulation results of dimensions (Figure S3) and noise
setups (Figures S4-S7).

Figure S3: Plots of relative error of $A$ relative to $A$ and $mA$ (dashed
blue curve), and relative error of $mA$ relative to $A$ and $mA$ (solid blue
curves). Plots of relative error of $A$ relative to $A$ and $hmA$ (dashed red
curves), and relative error of $hmA$ relative to $A$ and $hmA$ (solid red
curves).

Figure S4: Plots of winning proportion of $mA$ against $A$ where $A=$ PSO, BAT
and CSO categorized according to the dimension of the test functions. Each row
shows the result for one value of sigma $\sigma$, where
$\sigma=0.005,0.01,0.02,0.05$.

Figure S5: Plots of winning proportion of $hmA$ against $A$ where $A=$ PSO,
BAT and CSO categorized according to the dimension of the test functions. Each
row shows the result for one value of sigma $\sigma$, where
$\sigma=0.005,0.01,0.02,0.05$.

Figure S6: Plots of winning proportion of $mA-t$ against $A$ where $A=$ PSO,
BAT and CSO categorized according to the dimension of the test functions. Each
row shows the result for one value of degrees of freedom $df$, where
$df=0.005,0.01,0.02,0.05$.

Figure S7: Plots of winning proportion of $hmA-t$ against $A$ where $A=$ PSO,
BAT and CSO categorized according to the dimension of the test functions. Each
row shows the result for one value of degrees of freedom $df$, where
$df=0.005,0.01,0.02,0.05$.

##  Appendix C Differential Evolution (DE) and its modification

Here we review the DE algorithm and call the particles in DE as agents. In
each DE iteration, each agent considers a new proposed location and move to
the new location if it provides a better function value. The proposed movement
is obtained by truncating the difference between two other random agents at a
random component. More specifically, DE repeats steps 1 to 6 until the user-
specified termination criterion is met.

  1. 1.

Repeat the following step for each agent $\mathbf{x}_{i}(t)$

  2. 2.

Randomly find two other agents $\mathbf{x}_{j}(t)$ and $\mathbf{x}_{k}(t)$.

  3. 3.

Compute
$\mathbf{y}_{i}(t)=\mathbf{x}_{i}(t)+F(\mathbf{x}_{j}(t)-\mathbf{x}_{k}(t))$,
where $F\in[0,2]$ is a constant known as _differential weight_ and its default
value is $0.8$.

  4. 4.

Randomly find a component index $\ell\in\\{1,\ldots,d\\}$.

  5. 5.

For component $m\in\\{1,\ldots,d\\}\setminus\\{\ell\\}$, replace the $m$th
component of $\mathbf{y}_{i}(t)$ by the corresponding component of
$\mathbf{x}_{i}(t)$ with probability $1-\omega$, where $\omega\in[0,1]$ is a
constant known as _crossover probability_ , and its default value is 0.9.

  6. 6.

If $f(\mathbf{x}_{i}(t))>f(\mathbf{y}_{i}(t))$, let
$\mathbf{x}_{i}(t+1)=\mathbf{y}_{i}(t)$; otherwise let
$\mathbf{x}_{i}(t+1)=\mathbf{x}_{i}(t)$.

Our modified DE, denoted by mDE, is to add an additional step for step 6
6’. Let
$\mathbf{y}^{\prime}_{i}(t)=\chi_{\mathcal{D}}(\chi_{\mathcal{D}}(\mathbf{y}_{i}(t))+\mathbf{w}_{\ell}(t))$
and

$\mathbf{x}_{i}(t+1)=\begin{cases}\mathbf{y}^{\prime}_{i}(t),\quad\text{if
}f(\mathbf{y}^{\prime}_{i}(t))<f(\mathbf{x}_{i}(t));\\\
\mathbf{x}_{i}(t),\quad\text{if }f(\mathbf{y}^{\prime}_{i}(t))\geq
f(\mathbf{x}_{i}(t)).\\\ \end{cases}$ (6)

The simulations use the same setups as the ones for PSO. The results are
plotted in Figures S8 and S9. The modification actually worsen the
optimization output.

Figure S8: Plots of winning proportion of $mDE$ against $DE$ (left), and that
of $hmDE$ against $DE$ (right) where DE categorized according to the dimension
of the test functions.

Figure S9: Plots of relative error of $DE$ relative to $DE$ and $mDE$ (dashed
blue curve), and relative error of $mDE$ relative to $DE$ and $mDE$ (solid
blue curve). Plots of relative error of $DE$ relative to $DE$ and $hmDE$
(dashed red curves), and relative error of $hmDE$ relative to $DE$ and $hmDE$
(solid red curve).