9 Population Methods
Previous chapters have focused on methods where a single design point is moved
incrementally toward a minimum. This chapter presents a variety of population
methods that involve optimization using a collection of design points, called individuals. Having a large number of individuals distributed throughout the design
space can help the algorithm avoid becoming stuck in a local minimum. Information at different points in the design space can be shared between individuals to
globally optimize the objective function. Most population methods are stochastic
in nature, and it is generally easy to parallelize the computation.

9.1 Initialization
Population methods begin with an initial population, just as descent methods
require an initial design point. The initial population should be spread over thedesign space to increase the chances that the samples are close to the best regions.
This section presents several initialization methods, but more advanced sampling
methods are discussed in detail in chapter 13.
We can often constrain the design variables to a region of interest consisting ofa hyperrectangle defined by lower and upper bounds a and b. Initial populations
can be sampled from a uniform distribution for each coordinate:1
( j)

xi

∼ U ( a i , bi )(9.1)

where x( j) is the jth individual in the population as seen in algorithm 9.1.
Another common approach is to use a multivariate normal distribution centered over a region of interest. The covariance matrix is typically diagonal, withdiagonal entries scaled to adequately cover the search space. Algorithm 9.2 provides an implementation.

1
Some population methods require additional information to
be associated with the individual,such as velocity in the case of particle swarm optimization, discussed
later. Velocity is often initialized according to a uniform or normal distribution.

148

c ha p te r 9. p op u l a ti on me thod sfunction rand_population_uniform(m, a, b)
d = length(a)return [a+rand(d).*(b-a) for i in 1:m]
end

Algorithm 9.1. A method for sampling an initial population of m design points over a uniform hyperrectangle with lower-bound vector
a and upper-bound vector b.

using Distributions
function rand_population_normal(m, μ, Σ)
D = MvNormal(μ,Σ)
return [rand(D) for i in 1:m]
end

Algorithm 9.2. A method for sampling an initial population of m design points using a multivariate
normal distribution with mean μand covariance Σ.

Uniform and normal distributions limit the covered design space to a concentrated region. The Cauchy distribution (figure 9.1) has an unbounded varianceand can cover a much broader space. Algorithm 9.3 provides an implementation.
Figure 9.2, on the next page, compares example initial populations generatedusing different methods.
NormalCauchyp( x )

0.4

0.2

0

−5

0
x

5

9.2 Genetic Algorithms
Genetic algorithms (algorithm 9.4) borrow inspiration from biological evolution,
where fitter individuals are more likely to pass on their genes to the next generation.2 An individual’s fitness for reproduction is inversely related to the value ofthe objective function at that point. The design point associated with an individual
is represented as a chromosome. At each generation, the chromosomes of the fitter
individuals are passed on to the next generation after undergoing the genetic
operations of crossover and mutation.

Figure 9.1. A comparison of thenormal distribution with standard
deviation 1 and the Cauchy distribution with scale 1. Although
σ is sometimes used for the scaleparameter in the Cauchy distribution, this should not be confused with the standard deviation
since the standard deviation of the
Cauchy distribution is undefined.
The Cauchy distribution is heavytailed, allowing it to cover the design space more broadly.

D. E. Goldberg, Genetic Algorithms
in Search, Optimization, and Machine
Learning. Addison-Wesley, 1989.
2© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 2. g e n e ti c a l g ori thms

using Distributions
function rand_population_cauchy(m, μ, σ)n = length(μ)
return [[rand(Cauchy(μ[j],σ[j])) for j in 1:n] for i in 1:m]
end

Uniform42
x2

Algorithm 9.3. A method for sampling an initial population of m design points using a Cauchy distribution with location μ and scale
σ for each dimension. The location and scale are analogous to themean and standard deviation used
in a normal distribution.

Figure 9.2. Initial populations ofsize 1,000 sampled using a uniform
hyperrectangle with a = [−2, −2],b = [2, 2], a zero-mean normal distribution with diagonal covarianceΣ = I, and Cauchy distributions
centered at the origin with scale
σ = 1.

CauchyNormal

149

0

−2
−4
−4

−2

0
x1

2

4 −4

−2

0
x1

2

4 −4

−2

0
x1

2

4

function genetic_algorithm(f, population, k_max, S, C, M)
for k in 1 : k_max
parents = select(S, f.(population))children = [crossover(C,population[p[1]],population[p[2]])for p in parents]
population .= mutate.(Ref(M), children)end
population[argmin(f.(population))]end

Algorithm 9.4. The genetic
algorithm, which takes anobjective function f, an initial
population, number of iterations
k_max, a SelectionMethod S,
a CrossoverMethod C, and a
MutationMethod M.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

150

c ha p te r 9. p op u l a ti on me thod s

9.2.1 Chromosomes
There are several ways to represent chromosomes. The simplest is the binary
string chromosome, a representation that is similar to the way DNA is encoded.3 A
random binary string of length d can be generated using bitrand(d). A binary
string chromosome is depicted in figure 9.3.

Instead of a binary representation,
DNA contains four nucleobases:
adenine, thymine, cytosine, andguanine, which are often abbreviated A, T, C, and G.
3

Figure 9.3. A chromosome represented as a binary string.

Binary strings are often used due to the ease of expressing crossover andmutation. Unfortunately, the process of decoding a binary string and producing adesign point is not always straightforward. Sometimes the binary string might not
represent a valid point in the design space. It is often more natural to represent achromosome using a list of real values. Such real-valued chromosomes are vectors
in R d that directly correspond to points in the design space.

9.2.2 Initialization
Genetic algorithms start with a random initial population. Binary string chromosomes are typically initialized using random bit strings as seen in algorithm 9.5.
Real-valued chromosomes are typically initialized using the methods from theprevious section.
rand_population_binary(m, n) = [bitrand(n) for i in 1:m]

Algorithm 9.5. A method for sampling random starting populationsof m bit-string chromosomes oflength n.

9.2.3 Selection
Selection is the process of choosing chromosomes to use as parents for the next
generation. For a population with m chromosomes, a selection method will produce a list of m parental pairs4 for the m children of the next generation. The
selected pairs may contain duplicates.

In some cases, one might use
groups, should one wish to combine more than two parents to forma child.

4© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 2. g e n e ti c a l g ori thms

151

Figure 9.4. Truncation selection
with a population size m = 7 andsample size k = 3. The height of abar indicates its objective function
value whereas its color indicateswhat individual it corresponds to.

yyThere are several approaches for biasing the selection toward the fittest (algorithm 9.6). In truncation selection (figure 9.4), we sample parents from among thebest k chromosomes in the population. In tournament selection (figure 9.5), eachparent is the fittest out of k randomly chosen chromosomes of the population. In
roulette wheel selection (figure 9.6), also known as fitness proportionate selection, eachparent is chosen with a probability proportional to its performance relative to the
population. Since we are interested in minimizing an objective function f , the fitness of the ith individual x (i) is inversely related to y(i) = f ( x (i) ). There are different ways to transform a collection y(1) , . . . , y(m) into fitnesses. A simple approachis to assign the fitness of individual i according to max{y(1) , . . . , y(m) } − y(i) .

individualyyindividual

individualylikelihood

individual

individual

individual

Figure 9.5. Tournament selection
with a population size m = 7
and a sample size k = 3, which
is run separately for each parent.
The height of a bar indicates its objective function value whereas its
color indicates what individual itcorresponds to.

Figure 9.6. Roulette wheel selection with a population size m = 7,
which is run separately for eachparent. The approach used causes
the individual with the worst objective function value to have a zero
likelihood of being selected. The
height of a bar indicates its objective function value (left), or its likelihood (right), whereas its color
indicates what individual it corresponds to.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

152 c ha p te r 9. p op u l a ti on me thod s

abstract type SelectionMethod end# Pick pairs randomly from top k parents
struct TruncationSelection <: SelectionMethod
k # top k to keependfunction select(t::TruncationSelection, y)p = sortperm(y)return [p[rand(1:t.k, 2)] for i in y]
end

Algorithm 9.6. Several selection
methods for genetic algorithms. Calling select with a
SelectionMethod and a list of
objective function values y willproduce a list of parental pairs.

# Pick parents by choosing best among random subsets
struct TournamentSelection <: SelectionMethodkendfunction select(t::TournamentSelection, y)
getparent() = begin
p = randperm(length(y))
p[argmin(y[p[1:t.k]])]end
return [[getparent(), getparent()] for i in y]end# Sample parents proportionately to fitness
struct RouletteWheelSelection <: SelectionMethod end
function select(::RouletteWheelSelection, y)y = maximum(y) .- ycat = Categorical(normalize(y, 1))
return [rand(cat, 2) for i in y]end© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 2. g e n e ti c a l g ori thms

153

9.2.4 Crossover
Crossover combines the chromosomes of parents to form children. As with selection, there are several crossover schemes (algorithm 9.7).
• In single-point crossover (figure 9.7), the first portion of parent A’s chromosome
forms the first portion of the child chromosome, and the latter portion of parent
B’s chromosome forms the latter part of the child chromosome. The crossover
point where the transition occurs is determined uniformly at random.
parent A
parent B

Figure 9.7. Single-point crossover.

childcrossover point• In two-point crossover (figure 9.8), we use two random crossover points.
parent A
parent B
child

Figure 9.8. Two-point crossover.

crossover point 1

crossover point 2• In uniform crossover (figure 9.9), each bit has a fifty percent chance of comingfrom either one of the two parents. This scheme is equivalent to each point
having a fifty percent chance of being a crossover point.
parent A
parent B

Figure 9.9. Uniform crossover.

child

The previous crossover methods also work for real-valued chromosomes. We
can, however, define an additional crossover routine that interpolates between
real values (algorithm 9.8). Here, the real values are linearly interpolated between
the parents’ values xa and xb :x ← (1 − λ)xa + λxb(9.2)

where λ is a scalar parameter typically set to one-half.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

154 c ha p te r 9. p op u l a ti on me thod s

abstract type CrossoverMethod end
struct SinglePointCrossover <: CrossoverMethod end
function crossover(::SinglePointCrossover, a, b)
i = rand(1:length(a))
return vcat(a[1:i], b[i+1:end])
end
struct TwoPointCrossover <: CrossoverMethod end
function crossover(::TwoPointCrossover, a, b)
n = length(a)i, j = rand(1:n, 2)
if i > j(i,j) = (j,i)
endreturn vcat(a[1:i], b[i+1:j], a[j+1:n])end

Algorithm 9.7. Several crossover
methods for genetic algorithms. Calling crossover witha CrossoverMethod and twoparents a and b will produce achild chromosome that contains
a mixture of the parents’ genetic
codes. These methods work forboth binary string and real-valued
chromosomes.

struct UniformCrossover <: CrossoverMethod end
function crossover(::UniformCrossover, a, b)
child = copy(a)for i in 1 : length(a)if rand() < 0.5
child[i] = b[i]endendreturn childend

struct InterpolationCrossover <: CrossoverMethodλendcrossover(C::InterpolationCrossover, a, b) = (1-C.λ)*a + C.λ*b

Algorithm 9.8. A crossover
method for real-valued chromosomes which performs linear
interpolation between the parents.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 2. g e n e ti c a l g ori thms 155

9.2.5 Mutation
If new chromosomes were produced only through crossover, many traits that were
not present in the initial random population could never occur, and the most-fit
genes could saturate the population. Mutation allows new traits to spontaneously
appear, allowing the genetic algorithm to explore more of the state space. Childchromosomes undergo mutation after crossover.
Each bit in a binary-valued chromosome typically has a small probability of
being flipped (figure 9.10). For a chromosome with m bits, this mutation rate is
typically set to 1/m, yielding an average of one mutation per child chromosome.
Mutation for real-valued chromosomes can be implemented using bitwise flips,
but it is more common to add zero-mean Gaussian noise. Algorithm 9.9 provides
implementations.
before mutation
after mutation

Figure 9.10. Mutation for binary
string chromosomes gives each bita small probability of flipping.

abstract type MutationMethod end
struct BitwiseMutation <: MutationMethodλendfunction mutate(M::BitwiseMutation, child)return [rand() < M.λ ? !v : v for v in child]end

Algorithm 9.9. The bitwise mutation method for binary string chromosomes and the Gaussian mutation method for real-valued chromosomes. Here, λ is the mutation
rate, and σ is the standard deviation.

struct GaussianMutation <: MutationMethodσendfunction mutate(M::GaussianMutation, child)return child + randn(length(child))*M.σ
end

Figure 9.11 illustrates several generations of a genetic algorithm. Example 9.1
shows how to combine selection, crossover, and mutation strategies discussed inthis section.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

c ha p te r 9. p op u l a ti on me thod s

x2

156

x1

x1

x1

x1
Figure 9.11. A genetic algorithmwith truncation selection, single
point crossover, and Gaussian mutation with σ = 0.1 applied to theMichalewicz function defined inappendix B.5.

We will demonstrate using genetic algorithms to optimize a simple function.
import Random: seed!
import LinearAlgebra: norm
seed!(0) # set random seed for reproducible results
f = x->norm(x)m = 100 # population size
k_max = 10 # number of iterations
population = rand_population_uniform(m, [-3, 3], [3,3])
S = TruncationSelection(10) # select top 10
C = SinglePointCrossover()M = GaussianMutation(0.5) # small mutation rate
x = genetic_algorithm(f, population, k_max, S, C, M)@show xExample 9.1. Demonstration of using a genetic algorithm for optimizing a simple function.

x = [-0.00994906141228906, -0.05198433759424115]© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 3. d i ffe re n ti a l e vol u ti on

157

9.3 Differential Evolution
Differential evolution (algorithm 9.10) attempts to improve each individual in the
population by recombining other individuals in the population according to asimple formula.5 It is parameterized by a crossover probability p and a differential
weight w. Typically, w is between 0.4 and 1. For each individual x:
1. Choose three random distinct individuals a, b, and c.

S. Das and P. N. Suganthan, “Differential Evolution: A Survey of theState-of-the-Art,” IEEE Transactions
on Evolutionary Computation, vol. 15,
no. 1, pp. 4–31, 2011.
5

2. Construct an interim design z = a + w · (b − c) as shown in figure 9.12.

3. Choose a random dimension j ∈ [1, . . . , n] for optimization in n dimensions.

4. Construct the candidate individual x′ using binary crossover.
z if i = j or with probability pixi′ = xi otherwise

z

(9.3)

5. Insert the better design between x and x′ into the next generation.
The algorithm is demonstrated in figure 9.13.
using StatsBase
function differential_evolution(f, population, k_max; p=0.5, w=1)
n, m = length(population[1]), length(population)
for k in 1 : k_max
for (k,x) in enumerate(population)a, b, c = sample(population,
Weights([j!=k for j in 1:m]), 3, replace=false)z = a + w*(b-c)
j = rand(1:n)
x′ = [i == j || rand() < p ? z[i] : x[i] for i in 1:n]
if f(x′) < f(x)x[:] = x′
end
end
end
return population[argmin(f.(population))]endbw (b − c)ab−ccFigure 9.12. Differential evolution
takes three individuals a, b, and c
and combines them to form the candidate individual z.
Algorithm 9.10. Differential evolution, which takes an objective function f, a population population,
a number of iterations k_max, acrossover probability p, and a differential weight w. The best individual is returned.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

c ha p te r 9. p op u l a ti on me thod s

x2

x2

158

x1

x1

x1

9.4 Particle Swarm Optimization
Particle swarm optimization introduces momentum to accelerate convergence toward minima.6 Each individual, or particle, in the population keeps track of its
current position, velocity, and the best position it has seen so far (algorithm 9.11).
Momentum allows an individual to accumulate speed in a favorable direction,
independent of local perturbations.
mutable struct Particle
xv
x_bestend

x1
Figure 9.13. Differential evolution
with p = 0.5 and w = 0.2 applied
to Ackley’s function, defined in appendix B.1.
J. Kennedy, R. C. Eberhart, andY. Shi, Swarm Intelligence. Morgan
Kaufmann, 2001.
6

Algorithm 9.11. Each particle inparticle swarm optimization has aposition x and velocity v in design
space and keeps track of the best
design point found so far, x_best.

At each iteration, each individual is accelerated toward both the best position it
has seen and the best position found thus far by any individual. The acceleration
is weighted by a random term, with separate random numbers being generated© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 5 . fi re fl y a l g ori thm 159

for each acceleration. The update equations are:x (i ) ← x (i ) + v (i )

(9.4)


(i )





v(i) ← wv(i) + c1 r1 xbest − x (i) + c2 r2 xbest − x (i)



(9.5)

where xbest is the best location found so far over all particles; w, c1 , and c2 are parameters; and r1 and r2 are random numbers drawn from U (0, 1).7 Algorithm 9.12
provides an implementation. Figure 9.14 shows several iterations of the algorithm.

Algorithm 9.12. Particle swarm optimization, which takes an objective function f, a list of particles
population, a number of iterations
k_max, an inertia w, and momentum coefficients c1 and c2. The default values are those used by R.
Eberhart and J. Kennedy, “A New
Optimizer Using Particle Swarm
Theory,” in International Symposiumon Micro Machine and Human Science, 1995.

function particle_swarm_optimization(f, population, k_max;
w=1, c1=1, c2=1)
n = length(population[1].x)
x_best, y_best = copy(population[1].x_best), Inf
for P in populationy = f(P.x)
if y < y_best; x_best[:], y_best = P.x, y; end
end
for k in 1 : k_max
for P in population
r1, r2 = rand(n), rand(n)P.x += P.v
P.v = w*P.v + c1*r1.*(P.x_best - P.x) +c2*r2.*(x_best - P.x)
y = f(P.x)
if y < y_best; x_best[:], y_best = P.x, y; end
if y < f(P.x_best); P.x_best[:] = P.x; end
end
end
return population
end

9.5

A common strategy is to allow theinertia w to decay over time.
7

X.-S. Yang, Nature-Inspired Metaheuristic Algorithms. Luniver Press,
2008. Interestingly, male fireflies
flash to attract members of the
opposite sex, but females sometimes flash to attract males of otherspecies, which they then eat.
8

Firefly Algorithm

The firefly algorithm (algorithm 9.13) was inspired by the manner in which fireflies
flash their lights to attract mates.8 In the firefly algorithm, each individual in the
population is a firefly and can flash to attract other fireflies. At each iteration, allfireflies are moved toward all more attractive fireflies. A firefly a is moved toward
a firefly b with greater attraction according toa ← a + βI (kb − ak)(b − a) + αǫ

(9.6)

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

160 c ha p te r 9. p op u l a ti on me thod s

where I is the intensity of the attraction and β is the source intensity. A random
walk component is included as well, where ǫ is drawn from a zero-mean, unit
covariance multivariate Gaussian, and α scales the step size. The resulting update
is a random walk biased toward brighter fireflies.9
The intensity I decreases as the distance r between the two fireflies increases
and is defined to be 1 when r = 0. One approach is to model the intensity as apoint source radiating into space, in which case the intensity decreases according
to the inverse square law
1
I (r ) = 2(9.7)
rAlternatively, if a source is suspended in a medium that absorbs light, the intensity
will decrease according to an exponential decay
I (r ) = e−γr2Yang recommends β = 1 and α ∈[0, 1]. If β = 0, the behavior is arandom walk.

9

(9.8)

where γ is a the light absorption coefficient.10
We generally want to avoid equation (9.7) in practice due to the singularity at r = 0. A combination of the inverse square law and absorption can beapproximated with a Gaussian brightness drop-off:
I (r ) = e−γr

Figure 9.14. The particle swarm
method with w = 0.1, c1 = 0.25,
and c2 = 2 applied to Wheeler’s
Ridge, appendix B.7.

The distance between fireflies
ceases to matter as γ approaches
zero.

10

(9.9)

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 6 . c u c koo se a rc h 161

A firefly’s attraction is proportional to its performance. Attraction affects only
whether one fly is attracted to another fly, whereas intensity affects how muchthe less attractive fly moves. Figure 9.15 shows a few iterations of the algorithm.
using Distributions
function firefly(f, population, k_max;
β=1, α=0.1, brightness=r->exp(-r^2))

end

9.6

m = length(population[1])N = MvNormal(Matrix(1.0I, m, m))
for k in 1 : k_max
for a in population, b in population
if f(b) < f(a)r = norm(b-a)
a[:] += β*brightness(r)*(b-a) + α*rand(N)
end
end
end
return population[argmin([f(x) for x in population])]

Algorithm 9.13. The firefly algorithm, which takes an objectivefunction f, a population flies consisting of design points, a numberof iterations k_max, a source intensity β, a random walk step size α,
and an intensity function I. The
best design point is returned.

Figure 9.15. Firefly search with
α = 0.5, β = 1, and γ = 0.1 applied to the Branin function (appendix B.3).

Cuckoo Search

Cuckoo search (algorithm 9.14) is another nature-inspired algorithm named after
the cuckoo bird, which engages in a form of brood parasitism.11 Cuckoos lay theireggs in the nests of other birds, often birds of other species. When this occurs,
the host bird may detect the invasive egg and then destroy it or establish a new

X.-S. Yang and S. Deb, “Cuckoo
Search via Lévy Flights,” in World
Congress on Nature & Biologically Inspired Computing (NaBIC), 2009.
11© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

162 c ha p te r 9. p op u l a ti on me thod snest somewhere else. However, there is also a chance that the egg is accepted andraised by the host bird.12
In cuckoo search, each nest represents a design point. New design points can beproduced using Lévy flights from nests, which are random walks with step-lengths
from a heavy-tailed distribution. A new design point can replace a nest if it has abetter objective function value, which is analogous to cuckoo eggs replacing theeggs of birds of other species.
The core rules are:

Interestingly, an instinct of newly
hatched cuckoos is to knock othereggs or hatchlings (those belonging to the host bird) out of the nest.

12

1. A cuckoo will lay an egg in a randomly chosen nest.
2. The best nests with the best eggs will survive to the next generation.
3. Cuckoo eggs have a chance of being discovered by the host bird, in which casethe eggs are destroyed.
Cuckoo search relies on random flights to establish new nest locations. Theseflights start from an existing nest and then move randomly to a new location.
While we might be tempted to use a uniform or Gaussian distribution for the walk,these restrict the search to a relatively concentrated region. Instead, cuckoo searchuses a Cauchy distribution, which has a heavier tail. In addition, the Cauchy
distribution has been shown to be more representative of the movements of otheranimals in the wild.13 Figure 9.16 shows a few iterations of cuckoo search.
Other nature-inspired algorithms include the artificial bee colony, the gray wolf
optimizer, the bat algorithm, glowworm swarm optimization, intelligent water
drops, and harmony search.14 There has been some criticism of the proliferationof methods that make analogies to nature without fundamentally contributing
novel methods and understanding.15

9.7 Hybrid Methods
Many population methods perform well in global search, being able to avoid local
minima and finding the best regions of the design space. Unfortunately, thesemethods do not perform as well in local search in comparison to descent methods.
Several hybrid methods16 have been developed to extend population methods withdescent-based features to improve their performance in local search. There are
two general approaches to combining population methods with local searchtechniques:17

For example, a certain speciesof fruit fly explores its surroundings using Cauchy-like steps separated by 90° turns. A. M. Reynolds
and M. A. Frye, “Free-Flight Odor
Tracking in Drosophila is Consistent with an Optimal Intermittent Scale-Free Search,” PLoS ONE,vol. 2, no. 4, e354, 2007.
13

See, for example, D. Simon, Evolutionary Optimization Algorithms. Wiley, 2013.

14

This viewpoint is expressed byK. Sörensen, “Metaheuristics—the
Metaphor Exposed,” International
Transactions in Operational Research,
vol. 22, no. 1, pp. 3–18, 2015.
15

In the literature, these kinds oftechniques are also referred to as
memetic algorithms or genetic local
search.

16

K. W. C. Ku and M.-W. Mak, “Exploring the Effects of Lamarckian
and Baldwinian Learning in Evolving Recurrent Neural Networks,”
in IEEE Congress on Evolutionary
Computation(CEC), 1997.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

17

9. 7. hy bri d me thod susing Distributions
mutable struct Nest
x # positiony # value, f(x)endfunction cuckoo_search(f, population, k_max;
p_a=0.1, C=Cauchy(0,1))m, n = length(population), length(population[1].x)
a = round(Int, m*p_a)for k in 1 : k_max
i, j = rand(1:m), rand(1:m)
x = population[j].x + [rand(C) for k in 1 : n]y = f(x)if y < population[i].y
population[i].x[:] = x
population[i].y = y
end

163

Algorithm 9.14. Cuckoo search,which takes an objective functionf, an initial set of nests population,
a number of iterations k_max, percent of nests to abandon p_a, andflight distribution C. The flight distribution is typically a centered
Cauchy distribution.

p = sortperm(population, by=nest->nest.y, rev=true)
for i in 1 : a
j = rand(1:m-a)+a
population[p[i]] = Nest(population[p[j]].x +[rand(C) for k in 1 : n],f(population[p[i]].x)
)
end

end

end
return population

Figure 9.16. Cuckoo search applied to the Branin function (appendix B.3).
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

164

c ha p te r 9. p op u l a ti on me thod s• In Lamarckian learning, the population method is extended with a local searchmethod that locally improves each individual. The original individual and its
objective function value are replaced by the individual’s optimized counterpart
and its objective function value.
• In Baldwinian learning, the same local search method is applied to each individual, but the results are used only to update the individual’s perceived objectivefunction value. Individuals are not replaced but are merely associated with
optimized objective function values, which are not the same as their actual
objective function value. Baldwinian learning can help prevent premature
convergence.
The difference between these approaches is illustrated in example 9.2.
Consider optimizing f ( x ) = −e− x − 2e−( x−3) using a population of individuals initialized near x = 0.
22yLamarckian

Baldwinian

0

0

−1

−1

−2

−2

−2

0

2
x

Example 9.2. A comparison of theLamarckian and Baldwinian hybrid methods.

4

−2

0

2

4

x

A Lamarckian local search update applied to this population would movethe individuals toward the local minimum, reducing the chance that futureindividuals escape and find the global optimum near x = 3. A Baldwinian
approach will compute the same update but leaves the original designs
unchanged. The selection step will value each design according to its value
from a local search.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

9. 8 . su mma ry

165

9.8 Summary• Population methods use a collection of individuals in the design space to guideprogression toward an optimum.
• Genetic algorithms leverage selection, crossover, and mutations to produce
better subsequent generations.
• Differential evolution, particle swarm optimization, the firefly algorithm, andcuckoo search include rules and mechanisms for attracting design points tothe best individuals in the population while maintaining suitable state space
exploration.
• Population methods can be extended with local search approaches to improve
convergence.

9.9 Exercises
Exercise 9.1. What is the motivation behind the selection operation in genetic
algorithms?
Exercise 9.2. Why does mutation play such a fundamental role in genetic algorithms? How would we choose the mutation rate if we suspect there is a better
optimal solution?
Exercise 9.3. If we observe that particle swarm optimization results in fast convergence to a nonglobal minimum, how might we change the parameters of thealgorithm?

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com