>Business >Evolution strategies from the ground up in Python

### Evolution strategies from the ground up in Python

Evolution strategies is a stochastic global optimization algorithm.

It is an evolutionary algorithm connected to others, like the genetic algorithm, even though it is developed particularly for continuous function optimization.

In this guide, you will find out how to implement the evolution strategies optimization algorithm.

After going through this guide, you will be aware of:

• Evolution strategies being a stochastic global optimization algorithm that draws inspiration from the biological theory of evolution through natural selection.
• There is a standard terminology for Evolution Strategies and two typical variants of the algorithm referenced to as (mu, lambda)-ES and (mu + lambda)-ES
• How to implement and apply the Evolution strategies algorithm to continuous objective functions.

Tutorial Summarization

This tutorial is subdivided into three portions, which are:

1] Evolution strategies

2] Develop a (mu, lambda)-ES

3] Develop a (mu + lambda)-ES

Evolution Strategies

Evolution strategies, at times referenced to as Evolution Strategy (singular) or ES, is a stochastic global optimization algorithm.

The strategy was created in the 1960s to be implemented manually by engineers for minimal drag designs in a wind tunnel.

The family of algorithms referred to as Evolution Strategies (ES) were developed by Ingo Rechenberg and Hans-Paul Schwefel at the Technical University of Berlin in the mid 1960s.

Evolution strategies is a variant of evolutionary algorithms and draws inspiration from the biological theory of evolution through means of natural selection. Not like other evolutionary algorithms, it doesn’t leverage any form of crossover, rather, alteration of candidate solutions is restricted to mutation operators. In this fashion, Evolution Strategies might be thought of as a variant of parallel stochastic hill climbing.

The algorithm consists of a population of candidate solutions that are arbitrarily produced. Every iteration  of the algorithm consists of initially evaluating the population of solutions, then deleting all but a subset of the ideal solutions, referenced to as truncation selection. The remainder of the solutions (the parents) each are leveraged as the foundation for producing a number of new candidate solutions (mutation) that substitute or compete with the parents for a position in the population for consideration in the subsequent iteration of the algorithm (generation).

There are an array of variations of this process and a standard terminology to summarize the algorithm. The size of the population is referenced to as lambda and the number of parents chosen every iteration is referenced to as mu.

The number of children developed from every parent is calculated as (lambda / mu) and parameters should be selected so that the division has no remainder.

• mu: The number of parents chosen every iteration.
• Lambda: Size of the population
• Lambda/mu: Number of children produced from every chosen parent.

A bracket notation is leveraged to detail the algorithm configuration, example, (mu, lambda)-ES. For instance, if mu=5 and lambda=20, then it would be summarized as (5, 20)-ES. A comma (,) demarcating the mu and lambda parameters signifies that the children substitute the parents directly every iteration of the algorithm.

• (mu, lambda)-ES: A variant of evolution strategies where children substitute parents.

A plus (+) separation of the mu and lambda parameters signifies that the children and the parents together will define the population for the subsequent iteration.

• (mu + lambda)-ES: A variant of evolution strategies where children and parents are included to the population.

A stochastic hill climbing algorithm can be implemented as an Evolution Strategy and would possess the notation (1 + 1)-ES.

This is the similie or canonical ES algorithm and there are several extensions and variations detailed in the literature.

Now that we are acquainted with Evolution Strategies we can look into how to go about implementing the algorithm.

Develop a (mu, lambda)-ES

In this portion of the blog, we will produce a (mu, lambda)-ES, that is, a variant of the algorithm where children substitute parents.

To start with, let’s define a challenging problem as the foundation for implementation of the algorithm.

The Ackley function is not an instance of a multimodal objective function that has a singular global optima and several local optima in which a local search might get stuck.

As such, a global optimization strategy is needed. It is a two-dimensional objective function that possesses a global optima at [0,0] which evaluates to 0.0

The instance below implements the Ackley and develops a three-dimensional surface plot displaying the global optima and several local optima.

 123456789101112131415161718192021222324252627282930 # ackley multimodal functionfrom numpy import arangefrom numpy import expfrom numpy import sqrtfrom numpy import cosfrom numpy import efrom numpy import pifrom numpy import meshgridfrom matplotlib import pyplotfrom mpl_toolkits.mplot3d import Axes3D # objective functiondef objective(x, y):return -20.0 * exp(-0.2 * sqrt(0.5 * (x**2 + y**2))) – exp(0.5 * (cos(2 * pi * x) + cos(2 * pi * y))) + e + 20 # define range for inputr_min, r_max = -5.0, 5.0# sample input range uniformly at 0.1 incrementsxaxis = arange(r_min, r_max, 0.1)yaxis = arange(r_min, r_max, 0.1)# create a mesh from the axisx, y = meshgrid(xaxis, yaxis)# compute targetsresults = objective(x, y)# create a surface plot with the jet color schemefigure = pyplot.figure()axis = figure.gca(projection=’3d’)axis.plot_surface(x, y, results, cmap=’jet’)# show the plotpyplot.show()

Running the instance creates the surface plot of the Ackley function displaying the major number of local optima.

We will be producing arbitrary candidate solutions as well as altered variants of existing candidate solutions. It is critical that all candidate solutions are within the bounds of the search problem.

To accomplish this, we will produce a function to check if a candidate solution is within the bounds of the search and then discard it and produce another solution if it is not.

The in_bounds() function below will take a candidate solution (point) and the definition of the bounds of the search space (bounds) and return True if the solution is within the bounds of the search of false otherwise.

 12345678 # check if a point is within the bounds of the searchdef in_bounds(point, bounds):# enumerate all dimensions of the pointfor d in range(len(bounds)):# check if out of bounds for this dimensionif point[d] < bounds[d, 0] or point[d] > bounds[d, 1]:return Falsereturn True

We can then leverage this function when generating the initial population of ‘lam’ (example, lambda) random candidate solutions.

For instance:

 12345678 …# initial populationpopulation = list()for _ in range(lam):candidate = Nonewhile candidate is None or not in_bounds(candidate, bounds):candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0])population.append(candidate)

Then, we can iterate over a static number of iterations of the algorithm. Every iteration first consists of evaluating every candidate solution within the population.

We will calculate the scores and store them in a separate parallel list.

 123 …# evaluate fitness for the populationscores = [objective(c) for c in population]

Then, we are required to choose the “mu” parents with the best scores, lowest scores in this case, as we are minimizing the objective function.

We will perform this in two steps. To start with, we will rank the candidate solutions on the basis of their scores in ascending fashion so that the solution with the lowest score possesses a rank of 0, the next possesses a rank 1 and so on. We can leverage a double call of the argsort function to accomplish this.

We will then leverage the ranks and choose those parents that possess a rank below the value ‘mu’. This implies if mu is set to 5 to choose 5 parents, only those parents with a rank between 0 and 4 will be chosen.

 12345 …# rank scores in ascending orderranks = argsort(argsort(scores))# select the indexes for the top mu ranked solutionsselected = [i for i,_ in enumerate(ranks) if ranks[i] < mu]

We can then develop children for every chosen parent.

To start with, we must calculate the cumulative number of children to create per parent.

 123 …# calculate the number of children per parentn_children = int(lam / mu)

We can then iterate over each parent and develop altered variants of each.

We will develop children leveraging a similar technique leveraged in stochastic hill climbing. Particularly, every variable will be sampled leveraging a Gaussian distribution with the present value as the mean and the standard deviation furnished as a “step size” hyerparameter.

 123456 …# create children for parentfor _ in range(n_children):child = Nonewhile child is None or not in_bounds(child, bounds):child = population[i] + randn(len(bounds)) * step_size

We can also check if every chosen parent is better than the best solution observed thus far so that we can return the ideal solution at the conclusion of the search.

 12345 …# check if this parent is the best solution ever seenif scores[i] < best_eval:best, best_eval = population[i], scores[i]print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval))

The created children can be included to a list and we can substitute the population with the list of children at the conclusion of the algorithm iteration.

 123 …# replace population with childrenpopulation = children

We can connect all of this together, into a function referred to as es_comma() that carries out the comma variant of the Evolution strategy algorithm.

The function takes on the name of the objective function, the bounds with regards to the search space, the number of iterations, the step size, and the mu and lambda hyperparameters and returns the ideal solution identified during the search and its evaluation.

 123456789101112131415161718192021222324252627282930313233343536 # evolution strategy (mu, lambda) algorithmdef es_comma(objective, bounds, n_iter, step_size, mu, lam):best, best_eval = None, 1e+10# calculate the number of children per parentn_children = int(lam / mu)# initial populationpopulation = list()for _ in range(lam):candidate = Nonewhile candidate is None or not in_bounds(candidate, bounds):candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0])population.append(candidate)# perform the searchfor epoch in range(n_iter):# evaluate fitness for the populationscores = [objective(c) for c in population]# rank scores in ascending orderranks = argsort(argsort(scores))# select the indexes for the top mu ranked solutionsselected = [i for i,_ in enumerate(ranks) if ranks[i] < mu]# create children from parentschildren = list()for i in selected:# check if this parent is the best solution ever seenif scores[i] < best_eval:best, best_eval = population[i], scores[i]print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval))# create children for parentfor _ in range(n_children):child = Nonewhile child is None or not in_bounds(child, bounds):child = population[i] + randn(len(bounds)) * step_sizechildren.append(child)# replace population with childrenpopulation = childrenreturn [best, best_eval]

Then, we can go about applying this algorithm to our Ackley objective function.

We will execute the algorithm for 5,000 iterations and leverage a step size of 0.15 within the search space. We will leverage a population size (lambda) of 100 select 20 parents (mu). These hyperparameters were selected following a bit of trail and error.

At the conclusion of the search, we will report the best candidate solution identified during the search.

 1234567891011121314151617 …# seed the pseudorandom number generatorseed(1)# define range for inputbounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])# define the total iterationsn_iter = 5000# define the maximum step sizestep_size = 0.15# number of parents selectedmu = 20# the number of children generated by parentslam = 100# perform the evolution strategy (mu, lambda) searchbest, score = es_comma(objective, bounds, n_iter, step_size, mu, lam)print(‘Done!’)print(‘f(%s) = %f’ % (best, score))

Connecting this together, the complete instance of application of the comma version of the Evolution strategies algorithm to the Ackley objective function is detailed below.

 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980 # evolution strategy (mu, lambda) of the ackley objective functionfrom numpy import asarrayfrom numpy import expfrom numpy import sqrtfrom numpy import cosfrom numpy import efrom numpy import pifrom numpy import argsortfrom numpy.random import randnfrom numpy.random import randfrom numpy.random import seed # objective functiondef objective(v):x, y = vreturn -20.0 * exp(-0.2 * sqrt(0.5 * (x**2 + y**2))) – exp(0.5 * (cos(2 * pi * x) + cos(2 * pi * y))) + e + 20 # check if a point is within the bounds of the searchdef in_bounds(point, bounds):# enumerate all dimensions of the pointfor d in range(len(bounds)):# check if out of bounds for this dimensionif point[d] < bounds[d, 0] or point[d] > bounds[d, 1]:return Falsereturn True # evolution strategy (mu, lambda) algorithmdef es_comma(objective, bounds, n_iter, step_size, mu, lam):best, best_eval = None, 1e+10# calculate the number of children per parentn_children = int(lam / mu)# initial populationpopulation = list()for _ in range(lam):candidate = Nonewhile candidate is None or not in_bounds(candidate, bounds):candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0])population.append(candidate)# perform the searchfor epoch in range(n_iter):# evaluate fitness for the populationscores = [objective(c) for c in population]# rank scores in ascending orderranks = argsort(argsort(scores))# select the indexes for the top mu ranked solutionsselected = [i for i,_ in enumerate(ranks) if ranks[i] < mu]# create children from parentschildren = list()for i in selected:# check if this parent is the best solution ever seenif scores[i] < best_eval:best, best_eval = population[i], scores[i]print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval))# create children for parentfor _ in range(n_children):child = Nonewhile child is None or not in_bounds(child, bounds):child = population[i] + randn(len(bounds)) * step_sizechildren.append(child)# replace population with childrenpopulation = childrenreturn [best, best_eval]  # seed the pseudorandom number generatorseed(1)# define range for inputbounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])# define the total iterationsn_iter = 5000# define the maximum step sizestep_size = 0.15# number of parents selectedmu = 20# the number of children generated by parentslam = 100# perform the evolution strategy (mu, lambda) searchbest, score = es_comma(objective, bounds, n_iter, step_size, mu, lam)print(‘Done!’)print(‘f(%s) = %f’ % (best, score))

Running the instance reports the candidate solution and scores every time an improved solution is identified, then reports the ideal solution identified at the conclusion of the search.

Your outcomes might demonstrate variance provided the stochastic nature of the algorithm or evaluation procedure, or variations in numerical accuracy. Take up running the instance a few times and contrast the average outcome.

In this scenario, we can observe that about 22 enhancements to performance were observed during the search and the ideal solution is near to the optima.

Without a doubt, this solution can be furnished as a beginning point to a local search algorithm to be further refined, a usual practice when leveraging a global optimization algorithm such as ES.

 123456789101112131415161718192021222324 0, Best: f([-0.82977995 2.20324493]) = 6.912490, Best: f([-1.03232526 0.38816734]) = 4.492401, Best: f([-1.02971385 0.21986453]) = 3.689542, Best: f([-0.98361735 0.19391181]) = 3.407962, Best: f([-0.98189724 0.17665892]) = 3.297472, Best: f([-0.07254927 0.67931431]) = 3.296413, Best: f([-0.78716147 0.02066442]) = 2.982793, Best: f([-1.01026218 -0.03265665]) = 2.695163, Best: f([-0.08851828 0.26066485]) = 2.003254, Best: f([-0.23270782 0.04191618]) = 1.665184, Best: f([-0.01436704 0.03653578]) = 0.151617, Best: f([0.01247004 0.01582657]) = 0.067779, Best: f([0.00368129 0.00889718]) = 0.0297025, Best: f([ 0.00666975 -0.0045051 ]) = 0.0244933, Best: f([-0.00072633 -0.00169092]) = 0.00530211, Best: f([2.05200123e-05 1.51343187e-03]) = 0.00434315, Best: f([ 0.00113528 -0.00096415]) = 0.00427418, Best: f([ 0.00113735 -0.00030554]) = 0.00337491, Best: f([ 0.00048582 -0.00059587]) = 0.00219704, Best: f([-6.91643854e-04 -4.51583644e-05]) = 0.001971504, Best: f([ 2.83063223e-05 -4.60893027e-04]) = 0.001313725, Best: f([ 0.00032757 -0.00023643]) = 0.00115Done!f([ 0.00032757 -0.00023643]) = 0.001147

Now that we are acquainted with how to implement the comma variant of evolution strategies, let’s look at how we could go about implementing the plus version.

Develop a (mu + lambda)-ES

The plus variant of the Evolution Strategies algorithm is very much like the comma version.

The primary variation is that children and the parents comprise the population at the conclusion rather than only the children. This facilitates the parents to compete with the children for selection in the subsequent iteration of the algorithm.

This can have the outcome of increasingly greedy behaviour by the search algorithm and possibly premature convergence to local optima (suboptimal solutions). The advantage is that the algorithm is able to exploit good candidate solutions that were identified and concentrate intently on candidate solutions in this region, possibly identifying further improvements.

We can go about implementing the plus version of the algorithm through alteration of the function to include parents to the population when developing the children.

 123 …# keep the parentchildren.append(population[i])

The updated variant of the function with this addition, and with a new name es_plus(), is detailed here.

 1234567891011121314151617181920212223242526272829303132333435363738 # evolution strategy (mu + lambda) algorithmdef es_plus(objective, bounds, n_iter, step_size, mu, lam):best, best_eval = None, 1e+10# calculate the number of children per parentn_children = int(lam / mu)# initial populationpopulation = list()for _ in range(lam):candidate = Nonewhile candidate is None or not in_bounds(candidate, bounds):candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0])population.append(candidate)# perform the searchfor epoch in range(n_iter):# evaluate fitness for the populationscores = [objective(c) for c in population]# rank scores in ascending orderranks = argsort(argsort(scores))# select the indexes for the top mu ranked solutionsselected = [i for i,_ in enumerate(ranks) if ranks[i] < mu]# create children from parentschildren = list()for i in selected:# check if this parent is the best solution ever seenif scores[i] < best_eval:best, best_eval = population[i], scores[i]print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval))# keep the parentchildren.append(population[i])# create children for parentfor _ in range(n_children):child = Nonewhile child is None or not in_bounds(child, bounds):child = population[i] + randn(len(bounds)) * step_sizechildren.append(child)# replace population with childrenpopulation = childrenreturn [best, best_eval]

We can go about applying this version of the algorithm to the Ackley objective function with the same hyperparameters leveraged in the prior section.

The complete instance is detailed below:

 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081 # evolution strategy (mu + lambda) of the ackley objective functionfrom numpy import asarrayfrom numpy import expfrom numpy import sqrtfrom numpy import cosfrom numpy import efrom numpy import pifrom numpy import argsortfrom numpy.random import randnfrom numpy.random import randfrom numpy.random import seed # objective functiondef objective(v):x, y = vreturn -20.0 * exp(-0.2 * sqrt(0.5 * (x**2 + y**2))) – exp(0.5 * (cos(2 * pi * x) + cos(2 * pi * y))) + e + 20 # check if a point is within the bounds of the searchdef in_bounds(point, bounds):# enumerate all dimensions of the pointfor d in range(len(bounds)):# check if out of bounds for this dimensionif point[d] < bounds[d, 0] or point[d] > bounds[d, 1]:return Falsereturn True # evolution strategy (mu + lambda) algorithmdef es_plus(objective, bounds, n_iter, step_size, mu, lam):best, best_eval = None, 1e+10# calculate the number of children per parentn_children = int(lam / mu)# initial populationpopulation = list()for _ in range(lam):candidate = Nonewhile candidate is None or not in_bounds(candidate, bounds):candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0])population.append(candidate)# perform the searchfor epoch in range(n_iter):# evaluate fitness for the populationscores = [objective(c) for c in population]# rank scores in ascending orderranks = argsort(argsort(scores))# select the indexes for the top mu ranked solutionsselected = [i for i,_ in enumerate(ranks) if ranks[i] < mu]# create children from parentschildren = list()for i in selected:# check if this parent is the best solution ever seenif scores[i] < best_eval:best, best_eval = population[i], scores[i]print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval))# keep the parentchildren.append(population[i])# create children for parentfor _ in range(n_children):child = Nonewhile child is None or not in_bounds(child, bounds):child = population[i] + randn(len(bounds)) * step_sizechildren.append(child)# replace population with childrenpopulation = childrenreturn [best, best_eval] # seed the pseudorandom number generatorseed(1)# define range for inputbounds = asarray([[-5.0, 5.0], [-5.0, 5.0]])# define the total iterationsn_iter = 5000# define the maximum step sizestep_size = 0.15# number of parents selectedmu = 20# the number of children generated by parentslam = 100# perform the evolution strategy (mu + lambda) searchbest, score = es_plus(objective, bounds, n_iter, step_size, mu, lam)print(‘Done!’)print(‘f(%s) = %f’ % (best, score))

Running the instance reports the candidate solution and scores every time an improved solution is identified, then reports the ideal solution identified at the conclusion of the search.

Your outcomes may demonstrate variance provided the stochastic nature of the algorithm or evaluation procedure, or variations in numerical accuracy. Take up running the instance a few times and contrast the average outcome.

In this scenario, we can observe that about 24 enhancements to performance were observed during the search. We can also observe that an improved final solution was identified with an evaluation of 0.000532, contrasted to 0.001147 identified with the comma version on this objective function.

 1234567891011121314151617181920212223242526 0, Best: f([-0.82977995 2.20324493]) = 6.912490, Best: f([-1.03232526 0.38816734]) = 4.492401, Best: f([-1.02971385 0.21986453]) = 3.689542, Best: f([-0.96315064 0.21176994]) = 3.489422, Best: f([-0.9524528 -0.19751564]) = 3.392662, Best: f([-1.02643442 0.14956346]) = 3.247842, Best: f([-0.90172166 0.15791013]) = 3.170902, Best: f([-0.15198636 0.42080645]) = 3.084313, Best: f([-0.76669476 0.03852254]) = 3.063653, Best: f([-0.98979547 -0.01479852]) = 2.621383, Best: f([-0.10194792 0.33439734]) = 2.523533, Best: f([0.12633886 0.27504489]) = 2.243444, Best: f([-0.01096566 0.22380389]) = 1.554764, Best: f([0.16241469 0.12513091]) = 1.440685, Best: f([-0.0047592 0.13164993]) = 0.775115, Best: f([ 0.07285478 -0.0019298 ]) = 0.341566, Best: f([-0.0323925 -0.06303525]) = 0.329516, Best: f([0.00901941 0.0031937 ]) = 0.0295032, Best: f([ 0.00275795 -0.00201658]) = 0.00997109, Best: f([-0.00204732 0.00059337]) = 0.00615195, Best: f([-0.00101671 0.00112202]) = 0.00434555, Best: f([ 0.00020392 -0.00044394]) = 0.001392804, Best: f([3.86555110e-04 6.42776651e-05]) = 0.001114357, Best: f([ 0.00013889 -0.0001261 ]) = 0.00053Done!f([ 0.00013889 -0.0001261 ]) = 0.000532

This portion of the blog furnishes additional resources on the subject if you are seeking to delve deeper.

Papers

Evolution Strategies: A Comprehensive Introduction, 2002

Books

Essentials of Mathematics, 2011

Algorithms for Optimization, 2019.

Computational Intelligence: An introduction, 2007.

Articles

Evolution strategy, Wikipedia

Evolution strategies, Scholarpedia

Conclusion

In this guide, you found out how to implement the evolution strategies optimization algorithm.

Particularly, you learned:

• Evolution strategies being a stochastic global optimization algorithm that draws inspiration by the biological theory of evolution through natural selection.
• There is a standard terminology for Evolution strategies and two usual variants of the algorithm referenced to as (mu, lambda)-ES and (mu + lambda)-ES
• How to implement and apply the Evolution strategies algorithm to continuous objective functions.