Self-Adaptive Evolution Strategy (SAES)

The ($\mu$/$\rho$ , $\lambda$)-SA-ES and ($\mu$/$\rho$ + $\lambda$)-SA-ES algorithms

Init pop

$\forall$ gen

$\quad$ $\forall$ child

$\quad\quad$ 1. select $\rho$ parents

$\quad\quad$ 2. recombination of selected parents (if $\rho > 1$)

$\quad\quad$ 3. mutation of $\sigma$ (individual strategy) : $\sigma \leftarrow \sigma ~ e^{\tau \mathcal{N}(0,1)}$

$\quad\quad$ 4. mutation of $\boldsymbol{x}$ (objective param) : $\boldsymbol{x} \leftarrow \boldsymbol{x} + \sigma ~ \mathcal{N}(0,1)$

$\quad\quad$ 5. eval $f(\boldsymbol{x})$

$\quad$ Select next gen individuals

A Python inplementation



In [ ]:

    
%matplotlib inline

import numpy as np
import pandas as pd
import math

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib import colors

from IPython.display import display



In [ ]:

    
# Set the objective function
#from ailib.optimize.functions import sphere as func
from ailib.optimize.functions import sphere2d as func
#from ailib.optimize.functions import additive_gaussian_noise as noise
from ailib.optimize.functions import multiplicative_gaussian_noise as noise
#from ailib.optimize.functions import additive_poisson_noise as noise

#func.noise = noise      # Comment this line to use a deterministic objective function

xmin = func.bounds[0]   # TODO
xmax = func.bounds[1]   # TODO



In [ ]:

    
# Set the objective function

from ailib.optimize.functions import sphere2d as func
#from ailib.optimize.functions import rosen2d as func
#from ailib.optimize.functions import himmelblau2d as func
#from ailib.optimize.functions import rastrigin2d as func
#from ailib.optimize.functions import easom2d as func
#from ailib.optimize.functions import crossintray2d as func   # <- ERROR
#from ailib.optimize.functions import holder2d as func        # <- ERROR

#from ailib.optimize.functions import additive_gaussian_noise as noise
#from ailib.optimize.functions import multiplicative_gaussian_noise as noise
#from ailib.optimize.functions import additive_poisson_noise as noise

func.translation_vector = np.ones(shape=func.ndim) * 3.

#func.noise = noise      # Comment this line to use a deterministic objective function

xmin = func.bounds[0]   # TODO
xmax = func.bounds[1]   # TODO

($\mu$/1+$\lambda$)-$\sigma$-Self-Adaptation-ES



In [ ]:

    
"""
This is a simple Python implementation of the (mu/1 +, lambda)-sigmaSA-ES
"""

mu = 3    # mu: the number of parents
lmb = 12  # lambda: the number of children
rho = 1   # rho: number of parents per child
selection_operator = '+'
isotropic_mutation = True        # TODO

assert selection_operator in (',', '+')

d = 2    # number of dimension of the solution space

num_gen = 20

tau = 1./math.sqrt(2.*d)          # self-adaptation learning rate

# Init the population ##########################

# "pop" array layout:
# - the first mu lines contain parents
# - the next lambda lines contain children
# - the first column contains the individual's strategy (sigma)
# - the last column contains the individual's assess (f(x))
# - the other columns contain the individual value (x)

sigma_col_label = "sigma"
y_col_label = "y"

pop = pd.DataFrame(np.full([mu+lmb, d+2], np.nan),
                   columns=[sigma_col_label] + ["x" + str(d) for d in range(d)] + [y_col_label])  # TODO: ajouter colonnes "set" ou "is_parent", "parents" (?) et "gen"

all_indices = slice(None, None)
parent_indices = slice(0, mu)
children_indices = slice(mu, None)

sigma_col = 0
x_cols = slice(1, -1)
y_col = -1

pop.iloc[parent_indices, sigma_col] = 1.                                    # init the parents strategy to 1.0
#pop.iloc[parent_indices, x_cols] = np.random.normal(0., 1., size=[mu, d])   # init the parents value
pop.iloc[parent_indices, x_cols] = np.random.uniform(low=-10., high=10., size=[mu, d])   # init the parents value
pop.iloc[parent_indices, y_col] = func(pop.iloc[parent_indices, x_cols].values.T)  # evaluate parents

#print("Init")
#display(pop)

cmap = cm.gnuplot2 # magma

fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(20, 6))
ax3.set_xlabel('gen')
ax3.set_ylabel('y')
#pop.plot.scatter(x="x0", y="x1", c=matplotlib.colors.rgb2hex(cmap(0.)), ax=ax1)
#pop.plot.scatter(x=0, y="y", c=matplotlib.colors.rgb2hex(cmap(0.)), loglog=True, ax=ax2)

for gen in range(num_gen):
    
    # Parent selection #############################
    
    if rho == 1:
        
        # Each child is made from one randomly selected parent
        selected_parent_indices = np.random.randint(mu, size=lmb)
    
    elif rho == mu:
    
        # Recombine all parents for each child
        raise NotImplemented() # TODO
    
    elif 1 < rho < mu:
    
        # Recombine rho randomly selected parents for each child
        raise NotImplemented() # TODO
    
    else:
    
        raise ValueError()
        
    #print("Parent selection")
    #display(selected_parent_indices)
    
    # Recombination ################################

    if rho == 1:
        
        # Each child is made from one randomly selected parent
        pop.iloc[children_indices] = pop.iloc[selected_parent_indices].values
        
    elif rho == mu:
        
        # Recombine all parents for each child
        raise NotImplemented() # TODO
        
    elif 1 < rho < mu:
        
        # Recombine rho randomly selected parents for each child
        raise NotImplemented() # TODO
        
    else:
        
        raise ValueError()
    
    pop.iloc[children_indices, y_col] = np.nan
    
    #print("Recombination")
    #display(pop)
    
    # Mutate children's sigma ######################
    
    pop.iloc[children_indices, sigma_col] *= np.exp(tau * np.random.normal(size=lmb))
    
    #print("Mutate children's sigma")
    #display(pop)

    # Mutate children's value ######################
    
    sigma_array = np.tile(pop.iloc[children_indices, sigma_col], [d,1]).T    # TODO: <- CHECK THIS !!!
    random_array = np.random.normal(size=[lmb,d])
    pop.iloc[children_indices, x_cols] += sigma_array * random_array

    #print("Mutate children's value")
    #display(pop)
    
    # Evaluate children ############################
    
    pop.iloc[children_indices, y_col] = func(pop.iloc[children_indices, x_cols].values.T)

    #print("Evaluate children")
    #display(pop)
    
    # SAVES SHOULD BE DONE HERE...
    color_str = matplotlib.colors.rgb2hex(cmap(float(gen) / num_gen))
    pop.plot.scatter(x="x0", y="x1", c=color_str, ax=ax1)
    pop.plot.scatter(x="sigma", y="y", c=color_str, loglog=True, ax=ax2)
    ax3.semilogy(np.full(shape=pop.y.shape, fill_value=gen), pop.y, '.', color=color_str)

    # Select the best individuals ##################
    
    if selection_operator == ',':
        pop.iloc[parent_indices] = np.nan
    
    pop = pop.sort_values(by=[y_col_label], na_position='last').reset_index(drop=True)
    
    pop.iloc[children_indices] = np.nan
    
    #print("Select the best individuals")
    #display(pop)

print("Result")
display(pop.iloc[parent_indices])

Some explanations about $\sigma$ and $\tau$



In [ ]:

    
d = np.arange(1, 30)
y = 1./np.sqrt(2. * d)

plt.plot(d, y)
plt.xlabel("d (number of dimensions)", fontsize=18)
plt.ylabel(r"$\tau$ (self-adaptation learning rate)", fontsize=18);



In [ ]:

    
tau



In [ ]:

    
d = 2.



In [ ]:

    
import random

sigma_list = [1.]

for i in range(1000):
    sigma_list.append(sigma_list[-1] * math.exp(tau * random.normalvariate(0., 1.)))         # mutate sigma

#sigma = sigma * exp(tau*randn)         # mutate sigma

plt.loglog(sigma_list);



In [ ]:

    
x = np.linspace(-4, 4, 100)
y1 = np.exp(1./math.sqrt(1.*d) * x)
y2 = np.exp(1./math.sqrt(2.*d) * x)
y3 = np.exp(1./math.sqrt(3.*d) * x)
y4 = np.exp(1./(2.*d) * x)

plt.plot(x, y1, label="tau1")
plt.plot(x, y2, label="tau2")
plt.plot(x, y3, label="tau3")
plt.plot(x, y4, label="tau4")

plt.legend();



In [ ]:

    
tau1 = 1./math.sqrt(1.*d)
tau2 = 1./math.sqrt(2.*d)
tau3 = 1./math.sqrt(3.*d)
tau4 = 1./(2.*d)

x1 = np.exp(tau1 * np.random.normal(size=[100000]))
x2 = np.exp(tau2 * np.random.normal(size=[100000]))
x3 = np.exp(tau3 * np.random.normal(size=[100000]))
x4 = np.exp(tau4 * np.random.normal(size=[100000]))

bins = np.linspace(0, 10, 100)
plt.hist(x1, bins=bins, alpha=0.5, label=r"$\exp\left(\frac{1}{\sqrt{d}} \mathcal{N}(0,1)\right)$", lw=2, histtype='step')
plt.hist(x2, bins=bins, alpha=0.5, label=r"$\exp\left(\frac{1}{\sqrt{2d}} \mathcal{N}(0,1)\right)$", lw=2, histtype='step')
plt.hist(x3, bins=bins, alpha=0.5, label=r"$\exp\left(\frac{1}{\sqrt{3d}} \mathcal{N}(0,1)\right)$", lw=2, histtype='step')
plt.hist(x4, bins=bins, alpha=0.5, label=r"$\exp\left(\frac{1}{2d} \mathcal{N}(0,1)\right)$", lw=2, histtype='step')

plt.xlim(-0.25, 7)

plt.axvline(1, color='k', linestyle='dotted')

plt.xlabel("mutation factor", fontsize=18)
plt.ylabel(r"counts", fontsize=18)

plt.legend(fontsize='x-large');

Other inplementations

AILib

Import required modules



In [ ]:

    
# Init matplotlib

%matplotlib inline

import matplotlib
matplotlib.rcParams['figure.figsize'] = (8, 8)



In [ ]:

    
import numpy as np
import time

from ailib.optimize import SAES



In [ ]:

    
# Plot functions
from ailib.optimize.utils import plot_contour_2d_solution_space
from ailib.optimize.utils import plot_2d_solution_space

from ailib.optimize.utils import array_list_to_array
from ailib.optimize.utils import plot_fx_wt_iteration_number
from ailib.optimize.utils import plot_err_wt_iteration_number
from ailib.optimize.utils import plot_err_wt_execution_time
from ailib.optimize.utils import plot_err_wt_num_feval

Define the objective function



In [ ]:

    
## Objective function: Rosenbrock function (Scipy's implementation)
#func = scipy.optimize.rosen



In [ ]:

    
# Set the objective function

from ailib.optimize.functions import sphere2d as func
#from ailib.optimize.functions import rosen2d as func
#from ailib.optimize.functions import himmelblau2d as func
#from ailib.optimize.functions import rastrigin2d as func
#from ailib.optimize.functions import easom2d as func
#from ailib.optimize.functions import crossintray2d as func   # <- ERROR
#from ailib.optimize.functions import holder2d as func        # <- ERROR

#from ailib.optimize.functions import additive_gaussian_noise as noise
#from ailib.optimize.functions import multiplicative_gaussian_noise as noise
#from ailib.optimize.functions import additive_poisson_noise as noise

func.translation_vector = np.ones(shape=func.ndim) * 3.

#func.noise = noise      # Comment this line to use a deterministic objective function

xmin = func.bounds[0]   # TODO
xmax = func.bounds[1]   # TODO



In [ ]:

    
#%%time

saes = SAES()

func.do_eval_logs = True
func.reset_eval_counters()
func.reset_eval_logs()

res = saes.minimize(func, init_pop_mean=0., init_pop_std=10., polt=True)

func.do_eval_logs = False

eval_x_array = np.array(func.eval_logs_dict['x']).T
eval_error_array = np.array(func.eval_logs_dict['fx']) - func(func.arg_min)



In [ ]:

    
res



In [ ]:

    
#eval_x_array



In [ ]:

    
plot_contour_2d_solution_space(func,
                               xmin=xmin,
                               xmax=xmax,
                               xstar=res,
                               xvisited=eval_x_array,
                               title="SAES");



In [ ]:

    
plot_err_wt_num_feval(eval_error_array, x_log=True, y_log=True)



In [ ]:

    
%%time

eval_error_array_list = []

NUM_RUNS = 100

for run_index in range(NUM_RUNS):
    saes = SAES()

    func.do_eval_logs = True
    func.reset_eval_counters()
    func.reset_eval_logs()

    res = saes.minimize(func, init_pop_mean=0., init_pop_std=1., lmb=6)

    func.do_eval_logs = False

    eval_error_array = np.array(func.eval_logs_dict['fx']) - func(func.arg_min)

    print("x* =", res)
    
    eval_error_array_list.append(eval_error_array);



In [ ]:

    
plot_err_wt_num_feval(array_list_to_array(eval_error_array_list), x_log=True, y_log=False, plot_option="mean")

Octave/Matlab

An external Octave (Matlab) implementation is available there: https://homepages.fhv.at/hgb/downloads/mu_mu_I_lambda-ES.oct.

Mathematica

An external Mathematica implementation is available there: https://homepages.fhv.at/hgb/downloads/mu_mu_I_lambda-ES.mat