object --+    
          |    
OOOptimizer --+
              |
             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. See also options scaling_of_variables. 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.fcts.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.fcts.rosen can be replaced by any properly defined objective function, see cma.fcts 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})
(4_w,8)-CMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=234)
>>>
>>> # optimize the ellipsoid function
>>> es.optimize(cma.fcts.elli, verb_disp=1)
Iterat #Fevals   function value     axis ratio  sigma   minstd maxstd min:sec
    1       8 2.093015112685775e+04 1.0e+00 9.27e-01  9e-01  9e-01 0:0.0
    2      16 4.964814235917688e+04 1.1e+00 9.54e-01  9e-01  1e+00 0:0.0
    3      24 2.876682459926845e+05 1.2e+00 1.02e+00  9e-01  1e+00 0:0.0
  100     800 6.809045875281943e-01 1.3e+02 1.41e-02  1e-04  1e-02 0:0.2
  200    1600 2.473662150861846e-10 8.0e+02 3.08e-05  1e-08  8e-06 0:0.5
  233    1864 2.766344961865341e-14 8.6e+02 7.99e-07  8e-11  7e-08 0:0.6
>>>
>>> cma.pprint(es.result())
(array([ -1.98546755e-09,  -1.10214235e-09,   6.43822409e-11,
        -1.68621326e-11]),
 4.5119610261406537e-16,
 1666,
 1672,
 209,
 array([ -9.13545269e-09,  -1.45520541e-09,  -6.47755631e-11,
        -1.00643523e-11]),
 array([  3.20258681e-08,   3.15614974e-09,   2.75282215e-10,
         3.27482983e-11]))
>>> assert es.result()[1] < 1e-9
>>> help(es.result)
Help on method result in module cma:

result(self) method of cma.CMAEvolutionStrategy instance

return (xbest, f(xbest), evaluations_xbest, evaluations, iterations, pheno(xmean), effective_stds)

The optimization loop can also be written explicitly.

>>> import cma
>>> es = cma.CMAEvolutionStrategy(4 * [1], 1)
>>> while not es.stop():
...    X = es.ask()
...    es.tell(X, [cma.fcts.elli(x) for x in X])
...    es.disp()
<output omitted>

achieving the same result as above.

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

>>> import cma
>>> import numpy as np
>>> es = cma.CMAEvolutionStrategy(10 * [0.2], 0.5, {'bounds': [0, np.inf]})
>>> 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.fcts.somenan(x, cma.fcts.elli) # might return np.NaN
...         X.append(x)
...         fit.append(curr_fit)
...     es.tell(X, fit)
...     es.logger.add()
...     es.disp()
<output omitted>
>>>
>>> assert es.result()[1] < 1e-9
>>> assert es.result()[2] < 9000  # by internal termination
>>> # es.logger.plot()  # will plot data
>>> # cma.show()  # display plot window

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

>>> es = cma.CMAEvolutionStrategy(5 * [3], 1,
...                 {"transformation": [lambda x: x**2+2, None]})
>>> es.optimize(cma.fcts.rosen)
<output omitted>
>>> assert cma.fcts.rosen(es.result()[0]) < 1e-6 + 5.530760944396627e+02
>>> assert es.result()[2] < 3300

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):

>>> import cma
>>> es = cma.CMAEvolutionStrategy(4 * [0.2], 0.5, {'verb_disp': 0})
>>> es.logger.disp_header()  # to understand the print of disp
Iterat Nfevals  function value    axis ratio maxstd   minstd
>>> while not es.stop():
...     X = es.ask()
...     es.tell(X, [cma.fcts.sphere(x) for x in X])
...     es.logger.add()  # log current iteration
...     es.logger.disp([-1])  # display info for last iteration
1      8 2.72769793021748e+03 1.0e+00 4.05e-01 3.99e-01
2     16 6.58755537926063e+03 1.1e+00 4.00e-01 3.39e-01
<output ommitted>
193   1544 3.15195320957214e-15 1.2e+03 3.70e-08 3.45e-11
>>> 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), output is not displayed:

>>> import cma, numpy as np
>>>
>>> # restart with increasing population size (IPOP)
>>> bestever = cma.BestSolution()
>>> for lam in 10 * 2**np.arange(8):  # 10, 20, 40, 80, ..., 10 * 2**7
...     es = cma.CMAEvolutionStrategy('6 - 8 * np.random.rand(9)',  # 9-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.fcts.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.pprint(es.best.__dict__)
...
...     bestever.update(es.best)
...
...     # show a plot
...     # logger.plot();
...     if bestever.f < 1e-8:  # global optimum was hit
...         break
<output omitted>
>>> assert es.result()[1] < 1e-8

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):

try:
    import multiprocessing as mp
    import cma
    es = cma.CMAEvolutionStrategy(22 * [0.0], 1.0, {'maxiter':10})
    pool = mp.Pool(es.popsize)
    while not es.stop():
        X = es.ask()
        f_values = pool.map_async(cma.felli, X).get()
        # use chunksize parameter as es.popsize/len(pool)?
        es.tell(X, f_values)
        es.disp()
        es.logger.add()
except ImportError:
    pass

The final example shows how to resume:

>>> import cma, pickle
>>>
>>> es = cma.CMAEvolutionStrategy(12 * [0.1],  # a new instance, 12-D
...                               0.5)         # initial std sigma0
>>> es.optimize(cma.fcts.rosen, iterations=100)
>>> pickle.dump(es, open('saved-cma-object.pkl', 'wb'))
>>> print('saved')
>>> del es  # let's start fresh
>>>
>>> es = pickle.load(open('saved-cma-object.pkl', 'rb'))
>>> print('resumed')
>>> es.optimize(cma.fcts.rosen, verb_disp=200)
>>> assert es.result()[2] < 15000
>>> cma.pprint(es.result())

Details

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

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. The update is applied after the default update and only before the covariance matrix is decomposed, which limits the additional computational burden to be at most a factor of three (typically smaller). A typical speed up factor (number of f-evaluations) is between 1.1 and two.

References: Jastrebski and Arnold, CEC 2006, Glasmachers et al, GECCO 2010.

Selective mirroring is implemented with option CMA_mirrors in the method get_mirror(). Only the method ask_and_eval() (used by fmin) will then sample selectively mirrored vectors. 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.

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

>>> import cma
>>> es = cma.CMAEvolutionStrategy([0,0,0,0], 0.3)
>>> while not es.stop() and es.best.f > 1e-6:  # my_desired_target_f_value
...     X = es.ask()  # get list of new solutions
...     fit = [cma.fcts.rosen(x) for x in X]  # call function rosen with each solution
...     es.tell(X, fit)  # feed values

>>> import cma
>>> es = cma.CMAEvolutionStrategy(cma.np.random.randn(3), 1)
>>> x = cma.np.random.randn(3)
>>> assert cma.Mh.vequals_approximately(es.mean - (x - es.mean), es.get_mirror(x, preserve_length=True))
>>> x = es.ask(1)[0]
>>> vals = (es.get_mirror(x) - es.mean) / (x - es.mean)
>>> assert cma.Mh.equals_approximately(sum(vals), len(vals) * vals[0])

>>> import cma
>>> x0, sigma0 = 8*[10], 1  # 8-D
>>> es = cma.CMAEvolutionStrategy(x0, sigma0)
>>> while not es.stop():
...     X, fit = es.ask_and_eval(cma.fcts.elli)  # handles NaN with resampling
...     es.tell(X, fit)  # pass on fitness values
...     es.disp(20) # print every 20-th iteration
>>> print('terminated on ' + str(es.stop()))
<output omitted>

while not es.stop():
    es.tell(*es.ask_and_eval(cma.fcts.elli))

import cma
func = cma.fcts.elli  # choose objective function
es = cma.CMAEvolutionStrategy(cma.np.random.rand(10), 1)
while not es.stop():
   X = es.ask()
   es.tell(X, [func(x) for x in X])
es.result()  # where the result can be found

return::

(xbest, f(xbest), evaluations_xbest, evaluations, iterations,
    pheno(xmean), effective_stds)

import cma

# prepare
(x0, sigma0) = ... # initial values from previous trial
X = ... # list of generated solutions from a previous trial
f = ... # respective list of f-values

# resume
es = cma.CMAEvolutionStrategy(x0, sigma0)
es.feedForResume(X, f)

# continue with func as objective function
while not es.stop():
   X = es.ask()
   es.tell(X, [func(x) for x in X])

>>> import cma, numpy
>>> es = cma.CMAEvolutionStrategy(numpy.ones(10), 1)
>>> xx = numpy.random.randn(2, 10)
>>> d = es.mahalanobis_norm(es.gp.geno(xx[0]-xx[1]))