cma.evolution_strategy.CMAEvolutionStrategy

Calling Sequences

• es = CMAEvolutionStrategy(x0, sigma0)

• es = CMAEvolutionStrategy(x0, sigma0, opts)

• es = CMAEvolutionStrategy(x0, sigma0).optimize(objective_fct)

• res = CMAEvolutionStrategy(x0, sigma0,
                          opts).optimize(objective_fct).result
  

Arguments

x0

initial solution, starting point. x0 is given as "phenotype" which means, if:

opts = {'transformation': [transform, inverse]}

is given and inverse is None, the initial mean is not consistent with x0 in that transform(mean) does not equal to x0 unless transform(mean) equals mean.

sigma0

initial standard deviation. The problem variables should have been scaled, such that a single standard deviation on all variables is useful and the optimum is expected to lie within about x0 +- 3*sigma0. Often one wants to check for solutions close to the initial point. This allows, for example, for an easier check of consistency of the objective function and its interfacing with the optimizer. In this case, a much smaller sigma0 is advisable.

opts

options, a dictionary with optional settings, see class CMAOptions.

Main interface / usage

The interface is inherited from the generic OOOptimizer class (see also there). An object instance is generated from:

es = cma.CMAEvolutionStrategy(8 * [0.5], 0.2)

The least verbose interface is via the optimize method:

es.optimize(objective_func)
res = es.result

More verbosely, the optimization is done using the methods stop, ask, and tell:

while not es.stop():
    solutions = es.ask()
    es.tell(solutions, [cma.ff.rosen(s) for s in solutions])
    es.disp()
es.result_pretty()

where ask delivers new candidate solutions and tell updates the optim instance by passing the respective function values (the objective function cma.ff.rosen can be replaced by any properly defined objective function, see cma.ff for more examples).

To change an option, for example a termination condition to continue the optimization, call:

es.opts.set({'tolfacupx': 1e4})

The class CMAEvolutionStrategy also provides:

(solutions, func_values) = es.ask_and_eval(objective_func)

and an entire optimization can also be written like:

while not es.stop():
    es.tell(*es.ask_and_eval(objective_func))

Besides for termination criteria, in CMA-ES only the ranks of the func_values are relevant.

Attributes and Properties

• inputargs: passed input arguments
• inopts: passed options
• opts: actually used options, some of them can be changed any time via opts.set, see class CMAOptions
• popsize: population size lambda, number of candidate solutions returned by ask ()
• logger: a CMADataLogger instance utilized by optimize

Examples

Super-short example, with output shown:

>>> import cma
>>> # construct an object instance in 4-D, sigma0=1:
>>> es = cma.CMAEvolutionStrategy(4 * [1], 1, {'seed':234})
...      # doctest: +ELLIPSIS
(4_w,8)-aCMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=234...)

and optimize the ellipsoid function

>>> es.optimize(cma.ff.elli, verb_disp=1)  # doctest: +ELLIPSIS
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      8 2.09...
>>> assert len(es.result) == 8, es.result
>>> assert es.result[1] < 1e-9, es.result

The optimization loop can also be written explicitly:

>>> es = cma.CMAEvolutionStrategy(4 * [1], 1)  # doctest: +ELLIPSIS
(4_w,8)-aCMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=...
>>> while not es.stop():
...    X = es.ask()
...    es.tell(X, [cma.ff.elli(x) for x in X])
...    es.disp()  # doctest: +ELLIPSIS
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1      8 ...

achieving the same result as above.

An example with lower bounds (at zero) and handling infeasible solutions:

>>> import numpy as np
>>> es = cma.CMAEvolutionStrategy(10 * [0.2], 0.5,
...         {'bounds': [0, np.inf]})  #doctest: +ELLIPSIS
(5_w,...
>>> while not es.stop():
...     fit, X = [], []
...     while len(X) < es.popsize:
...         curr_fit = None
...         while curr_fit in (None, np.nan):
...             x = es.ask(1)[0]
...             curr_fit = cma.ff.somenan(x, cma.ff.elli) # might return np.nan
...         X.append(x)
...         fit.append(curr_fit)
...     es.tell(X, fit)
...     es.logger.add()
...     es.disp()  #doctest: +ELLIPSIS
Itera...
>>>
>>> assert es.result[1] < 1e-9, es.result
>>> assert es.result[2] < 9000, es.result  # by internal termination
>>> # es.logger.plot()  # will plot data
>>> # cma.s.figshow()  # display plot window

An example with user-defined transformation, in this case to realize a lower bound of 2.

>>> import warnings
>>> with warnings.catch_warnings(record=True) as warns:
...     es = cma.CMAEvolutionStrategy(5 * [3], 0.1,
...                 {"transformation": [lambda x: x**2+1.2, None],
...                  "ftarget": 1e-7 + 5.54781521192,
...                  "verbose": -2,})
>>> warns[0].message  # doctest:+ELLIPSIS
UserWarning('in class GenoPheno: user defined transformations have not been tested thoroughly (...
>>> warns[1].message  # doctest:+ELLIPSIS
UserWarning('computed initial point...
>>> es.optimize(cma.ff.rosen, verb_disp=0)  #doctest: +ELLIPSIS
<cma...
>>> assert cma.ff.rosen(es.result[0]) < 1e-7 + 5.54781521192, es.result
>>> assert es.result[2] < 3300, es.result

The inverse transformation is (only) necessary if the BoundPenalty boundary handler is used at the same time.

The CMAEvolutionStrategy class also provides a default logger (cave: files are overwritten when the logger is used with the same filename prefix):

>>> es = cma.CMAEvolutionStrategy(4 * [0.2], 0.5, {'verb_disp': 0})
>>> es.logger.disp_header()  # annotate the print of disp
Iterat Nfevals  function value    axis ratio maxstd  minstd
>>> while not es.stop():
...     X = es.ask()
...     es.tell(X, [cma.ff.sphere(x) for x in X])
...     es.logger.add()  # log current iteration
...     es.logger.disp([-1])  # display info for last iteration   #doctest: +ELLIPSIS
    1  ...
>>> es.logger.disp_header()
Iterat Nfevals  function value    axis ratio maxstd  minstd
>>> # es.logger.plot() # will make a plot

Example implementing restarts with increasing popsize (IPOP):

>>> bestever = cma.optimization_tools.BestSolution()
>>> for lam in 10 * 2**np.arange(8):  # 10, 20, 40, 80, ..., 10 * 2**7
...     es = cma.CMAEvolutionStrategy(6 - 8 * np.random.rand(4),  # 4-D
...                                   5,  # initial std sigma0
...                                   {'popsize': lam,  # options
...                                    'verb_append': bestever.evalsall})
...     # logger = cma.CMADataLogger().register(es, append=bestever.evalsall)
...     while not es.stop():
...         X = es.ask()    # get list of new solutions
...         fit = [cma.ff.rastrigin(x) for x in X]  # evaluate each solution
...         es.tell(X, fit) # besides for termination only the ranking in fit is used
...
...         # display some output
...         # logger.add()  # add a "data point" to the log, writing in files
...         es.disp()  # uses option verb_disp with default 100
...
...     print('termination:', es.stop())
...     cma.s.pprint(es.best.__dict__)
...
...     bestever.update(es.best)
...
...     # show a plot
...     # logger.plot();
...     if bestever.f < 1e-8:  # global optimum was hit
...         break  #doctest: +ELLIPSIS
(5_w,...
>>> assert es.result[1] < 1e-8, es.result

On the Rastrigin function, usually after five restarts the global optimum is located.

Using the multiprocessing module, we can evaluate the function in parallel with a simple modification of the example (however multiprocessing seems not always reliable):

>>> from cma.fitness_functions import elli  # cannot be an instance method
>>> from cma.optimization_tools import EvalParallel2
>>> es = cma.CMAEvolutionStrategy(22 * [0.0], 1.0, {'maxiter':10})  # doctest:+ELLIPSIS
(6_w,13)-aCMA-ES (mu_w=...
>>> with EvalParallel2(elli, es.popsize + 1) as eval_all:
...     while not es.stop():
...         X = es.ask()
...         es.tell(X, eval_all(X))
...         es.disp()
...         # es.logger.add()  # doctest:+ELLIPSIS
Iterat...

The final example shows how to resume:

>>> import pickle
>>>
>>> es0 = cma.CMAEvolutionStrategy(12 * [0.1],  # a new instance, 12-D
...                                0.12)         # initial std sigma0
...   #doctest: +ELLIPSIS
(5_w,...
>>> es0.optimize(cma.ff.rosen, iterations=100)  #doctest: +ELLIPSIS
I...
>>> s = es0.pickle_dumps()  # return pickle.dumps(es) with safeguards
>>> # save string s to file like open(filename, 'wb').write(s)
>>> del es0  # let's start fresh
>>> # s = open(filename, 'rb').read()  # load string s from file
>>> es = pickle.loads(s)  # read back es instance from string
>>> # resuming
>>> es.optimize(cma.ff.rosen, verb_disp=200)  #doctest: +ELLIPSIS
  200 ...
>>> assert es.result[2] < 15000, es.result
>>> assert cma.s.Mh.vequals_approximately(es.result[0], 12 * [1], 1e-5), es.result
>>> assert len(es.result) == 8, es.result

Details

The following two enhancements are implemented, the latter is only turned on by default for very small population sizes.

Active CMA is implemented with option CMA_active and conducts an update of the covariance matrix with negative weights. The negative update is implemented, such that positive definiteness is guarantied. A typical speed up factor (number of f-evaluations) is between 1.1 and two.

References: Jastrebski and Arnold, Improving evolution strategies through active covariance matrix adaptation, CEC 2006. Hansen, The CMA evolution strategy: a tutorial, arXiv 2016.

Selective mirroring is implemented with option CMA_mirrors in the method get_mirror and get_selective_mirrors. The method ask_and_eval (used by fmin) will then sample selectively mirrored vectors within the iteration (CMA_mirrormethod==1). Otherwise, or if CMA_mirromethod==2, selective mirrors are injected for the next iteration. In selective mirroring, only the worst solutions are mirrored. With the default small number of mirrors, pairwise selection (where at most one of the two mirrors contribute to the update of the distribution mean) is implicitly guarantied under selective mirroring and therefore not explicitly implemented.

Update: pairwise selection for injected mirrors is also applied in the covariance matrix update: for all injected solutions, as for those from TPA, this is now implemented in that the recombination weights are constrained to be nonnegative for injected solutions in the covariance matrix (otherwise recombination weights are anyway nonnegative). This is a precaution to prevent failure when injected solutions are systematically bad (see e.g. https://github.com/CMA-ES/pycma/issues/124), but may not be "optimal" for mirrors.

References: Brockhoff et al, PPSN 2010, Auger et al, GECCO 2011.