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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | # ackley multimodal function from numpy import arange from numpy import exp from numpy import sqrt from numpy import cos from numpy import e from numpy import pi from numpy import meshgrid from matplotlib import pyplot from mpl_toolkits.mplot3d import Axes3D
# objective function def 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 input r_min, r_max = -5.0, 5.0 # sample input range uniformly at 0.1 increments xaxis = arange(r_min, r_max, 0.1) yaxis = arange(r_min, r_max, 0.1) # create a mesh from the axis x, y = meshgrid(xaxis, yaxis) # compute targets results = objective(x, y) # create a surface plot with the jet color scheme figure = pyplot.figure() axis = figure.gca(projection=’3d’) axis.plot_surface(x, y, results, cmap=’jet’) # show the plot pyplot.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.
| 1 2 3 4 5 6 7 8 | # check if a point is within the bounds of the search def in_bounds(point, bounds): # enumerate all dimensions of the point for d in range(len(bounds)): # check if out of bounds for this dimension if point[d] < bounds[d, 0] or point[d] > bounds[d, 1]: return False return True |
We can then leverage this function when generating the initial population of ‘lam’ (example, lambda) random candidate solutions.
For instance:
| 1 2 3 4 5 6 7 8 | … # initial population population = list() for _ in range(lam): candidate = None while 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.
| 1 2 3 | … # evaluate fitness for the population scores = [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.
| 1 2 3 4 5 | … # rank scores in ascending order ranks = argsort(argsort(scores)) # select the indexes for the top mu ranked solutions selected = [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.
| 1 2 3 | … # calculate the number of children per parent n_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.
| 1 2 3 4 5 6 | … # create children for parent for _ in range(n_children): child = None while 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.
| 1 2 3 4 5 | … # check if this parent is the best solution ever seen if 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.
| 1 2 3 | … # replace population with children population = 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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | # evolution strategy (mu, lambda) algorithm def es_comma(objective, bounds, n_iter, step_size, mu, lam): best, best_eval = None, 1e+10 # calculate the number of children per parent n_children = int(lam / mu) # initial population population = list() for _ in range(lam): candidate = None while candidate is None or not in_bounds(candidate, bounds): candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0]) population.append(candidate) # perform the search for epoch in range(n_iter): # evaluate fitness for the population scores = [objective(c) for c in population] # rank scores in ascending order ranks = argsort(argsort(scores)) # select the indexes for the top mu ranked solutions selected = [i for i,_ in enumerate(ranks) if ranks[i] < mu] # create children from parents children = list() for i in selected: # check if this parent is the best solution ever seen if scores[i] < best_eval: best, best_eval = population[i], scores[i] print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval)) # create children for parent for _ in range(n_children): child = None while child is None or not in_bounds(child, bounds): child = population[i] + randn(len(bounds)) * step_size children.append(child) # replace population with children population = children return [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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | … # seed the pseudorandom number generator seed(1) # define range for input bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]]) # define the total iterations n_iter = 5000 # define the maximum step size step_size = 0.15 # number of parents selected mu = 20 # the number of children generated by parents lam = 100 # perform the evolution strategy (mu, lambda) search best, 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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | # evolution strategy (mu, lambda) of the ackley objective function from numpy import asarray from numpy import exp from numpy import sqrt from numpy import cos from numpy import e from numpy import pi from numpy import argsort from numpy.random import randn from numpy.random import rand from numpy.random import seed
# objective function def objective(v): x, y = v 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
# check if a point is within the bounds of the search def in_bounds(point, bounds): # enumerate all dimensions of the point for d in range(len(bounds)): # check if out of bounds for this dimension if point[d] < bounds[d, 0] or point[d] > bounds[d, 1]: return False return True
# evolution strategy (mu, lambda) algorithm def es_comma(objective, bounds, n_iter, step_size, mu, lam): best, best_eval = None, 1e+10 # calculate the number of children per parent n_children = int(lam / mu) # initial population population = list() for _ in range(lam): candidate = None while candidate is None or not in_bounds(candidate, bounds): candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0]) population.append(candidate) # perform the search for epoch in range(n_iter): # evaluate fitness for the population scores = [objective(c) for c in population] # rank scores in ascending order ranks = argsort(argsort(scores)) # select the indexes for the top mu ranked solutions selected = [i for i,_ in enumerate(ranks) if ranks[i] < mu] # create children from parents children = list() for i in selected: # check if this parent is the best solution ever seen if scores[i] < best_eval: best, best_eval = population[i], scores[i] print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval)) # create children for parent for _ in range(n_children): child = None while child is None or not in_bounds(child, bounds): child = population[i] + randn(len(bounds)) * step_size children.append(child) # replace population with children population = children return [best, best_eval]
# seed the pseudorandom number generator seed(1) # define range for input bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]]) # define the total iterations n_iter = 5000 # define the maximum step size step_size = 0.15 # number of parents selected mu = 20 # the number of children generated by parents lam = 100 # perform the evolution strategy (mu, lambda) search best, 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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | 0, Best: f([-0.82977995 2.20324493]) = 6.91249 0, Best: f([-1.03232526 0.38816734]) = 4.49240 1, Best: f([-1.02971385 0.21986453]) = 3.68954 2, Best: f([-0.98361735 0.19391181]) = 3.40796 2, Best: f([-0.98189724 0.17665892]) = 3.29747 2, Best: f([-0.07254927 0.67931431]) = 3.29641 3, Best: f([-0.78716147 0.02066442]) = 2.98279 3, Best: f([-1.01026218 -0.03265665]) = 2.69516 3, Best: f([-0.08851828 0.26066485]) = 2.00325 4, Best: f([-0.23270782 0.04191618]) = 1.66518 4, Best: f([-0.01436704 0.03653578]) = 0.15161 7, Best: f([0.01247004 0.01582657]) = 0.06777 9, Best: f([0.00368129 0.00889718]) = 0.02970 25, Best: f([ 0.00666975 -0.0045051 ]) = 0.02449 33, Best: f([-0.00072633 -0.00169092]) = 0.00530 211, Best: f([2.05200123e-05 1.51343187e-03]) = 0.00434 315, Best: f([ 0.00113528 -0.00096415]) = 0.00427 418, Best: f([ 0.00113735 -0.00030554]) = 0.00337 491, Best: f([ 0.00048582 -0.00059587]) = 0.00219 704, Best: f([-6.91643854e-04 -4.51583644e-05]) = 0.00197 1504, Best: f([ 2.83063223e-05 -4.60893027e-04]) = 0.00131 3725, Best: f([ 0.00032757 -0.00023643]) = 0.00115 Done! 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.
| 1 2 3 | … # keep the parent children.append(population[i]) |
The updated variant of the function with this addition, and with a new name es_plus(), is detailed here.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 | # evolution strategy (mu + lambda) algorithm def es_plus(objective, bounds, n_iter, step_size, mu, lam): best, best_eval = None, 1e+10 # calculate the number of children per parent n_children = int(lam / mu) # initial population population = list() for _ in range(lam): candidate = None while candidate is None or not in_bounds(candidate, bounds): candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0]) population.append(candidate) # perform the search for epoch in range(n_iter): # evaluate fitness for the population scores = [objective(c) for c in population] # rank scores in ascending order ranks = argsort(argsort(scores)) # select the indexes for the top mu ranked solutions selected = [i for i,_ in enumerate(ranks) if ranks[i] < mu] # create children from parents children = list() for i in selected: # check if this parent is the best solution ever seen if scores[i] < best_eval: best, best_eval = population[i], scores[i] print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval)) # keep the parent children.append(population[i]) # create children for parent for _ in range(n_children): child = None while child is None or not in_bounds(child, bounds): child = population[i] + randn(len(bounds)) * step_size children.append(child) # replace population with children population = children return [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:
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 | # evolution strategy (mu + lambda) of the ackley objective function from numpy import asarray from numpy import exp from numpy import sqrt from numpy import cos from numpy import e from numpy import pi from numpy import argsort from numpy.random import randn from numpy.random import rand from numpy.random import seed
# objective function def objective(v): x, y = v 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
# check if a point is within the bounds of the search def in_bounds(point, bounds): # enumerate all dimensions of the point for d in range(len(bounds)): # check if out of bounds for this dimension if point[d] < bounds[d, 0] or point[d] > bounds[d, 1]: return False return True
# evolution strategy (mu + lambda) algorithm def es_plus(objective, bounds, n_iter, step_size, mu, lam): best, best_eval = None, 1e+10 # calculate the number of children per parent n_children = int(lam / mu) # initial population population = list() for _ in range(lam): candidate = None while candidate is None or not in_bounds(candidate, bounds): candidate = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] – bounds[:, 0]) population.append(candidate) # perform the search for epoch in range(n_iter): # evaluate fitness for the population scores = [objective(c) for c in population] # rank scores in ascending order ranks = argsort(argsort(scores)) # select the indexes for the top mu ranked solutions selected = [i for i,_ in enumerate(ranks) if ranks[i] < mu] # create children from parents children = list() for i in selected: # check if this parent is the best solution ever seen if scores[i] < best_eval: best, best_eval = population[i], scores[i] print(‘%d, Best: f(%s) = %.5f’ % (epoch, best, best_eval)) # keep the parent children.append(population[i]) # create children for parent for _ in range(n_children): child = None while child is None or not in_bounds(child, bounds): child = population[i] + randn(len(bounds)) * step_size children.append(child) # replace population with children population = children return [best, best_eval]
# seed the pseudorandom number generator seed(1) # define range for input bounds = asarray([[-5.0, 5.0], [-5.0, 5.0]]) # define the total iterations n_iter = 5000 # define the maximum step size step_size = 0.15 # number of parents selected mu = 20 # the number of children generated by parents lam = 100 # perform the evolution strategy (mu + lambda) search best, 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.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | 0, Best: f([-0.82977995 2.20324493]) = 6.91249 0, Best: f([-1.03232526 0.38816734]) = 4.49240 1, Best: f([-1.02971385 0.21986453]) = 3.68954 2, Best: f([-0.96315064 0.21176994]) = 3.48942 2, Best: f([-0.9524528 -0.19751564]) = 3.39266 2, Best: f([-1.02643442 0.14956346]) = 3.24784 2, Best: f([-0.90172166 0.15791013]) = 3.17090 2, Best: f([-0.15198636 0.42080645]) = 3.08431 3, Best: f([-0.76669476 0.03852254]) = 3.06365 3, Best: f([-0.98979547 -0.01479852]) = 2.62138 3, Best: f([-0.10194792 0.33439734]) = 2.52353 3, Best: f([0.12633886 0.27504489]) = 2.24344 4, Best: f([-0.01096566 0.22380389]) = 1.55476 4, Best: f([0.16241469 0.12513091]) = 1.44068 5, Best: f([-0.0047592 0.13164993]) = 0.77511 5, Best: f([ 0.07285478 -0.0019298 ]) = 0.34156 6, Best: f([-0.0323925 -0.06303525]) = 0.32951 6, Best: f([0.00901941 0.0031937 ]) = 0.02950 32, Best: f([ 0.00275795 -0.00201658]) = 0.00997 109, Best: f([-0.00204732 0.00059337]) = 0.00615 195, Best: f([-0.00101671 0.00112202]) = 0.00434 555, Best: f([ 0.00020392 -0.00044394]) = 0.00139 2804, Best: f([3.86555110e-04 6.42776651e-05]) = 0.00111 4357, Best: f([ 0.00013889 -0.0001261 ]) = 0.00053 Done! f([ 0.00013889 -0.0001261 ]) = 0.000532 |
Further Reading
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.
