CGO 09-1: Meta-heuristics

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits import mplot3d
import cma

CMA-ES Example

Small example of optimization of the Rastrigin function in 2-D. This uses the pycma library developed and maintained by the author of the CMA-ES algorithm Nikolaus Hansen.

Installation

pip install cma

Rastrigin Function

The Rastrigin function is characterized by the following equation for a vector $\mathbf{x} \in \mathbb{R}^n$:

\[f(\mathbf{x}) = A n + \sum_{i=1}^n \left(x_i^2 - 10\cos(2 \pi x_i)\right)\]

It is a very challenging function for global optimization methods as it has many local minimums. The function has a single global minimum located at 0.

def f(x): # Rastrigin function
    return (10 + x**2 - 10 * np.cos( 2*np.pi*x )).sum()

# Visualization of the Rastrigin function
RESOLUTION = 128
RANGE = 5.12
x = np.linspace( -RANGE, RANGE, RESOLUTION )
y = np.linspace( -RANGE, RANGE, RESOLUTION )
X,Y = np.meshgrid(x, y)
Z = np.zeros( (RESOLUTION, RESOLUTION) )
for i in range(Z.shape[0]):
    for j in range(Z.shape[1]):
        Z[i,j] = f( np.array( (X[i,j], Y[i,j]) ) )

# 2D Plot
plt.contourf(X, Y, Z, 20, cmap=cm.coolwarm)
plt.axis('equal')
plt.colorbar()

# 3D Plot
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, cmap=cm.coolwarm )
plt.show()

x0 = 2 * [5]
C0 = 10
es = cma.CMAEvolutionStrategy( x0, C0 )

# Using "Show and Tell" interface
while not es.stop():
    solutions = es.ask()
    es.tell(solutions, [f(x) for x in solutions])
    es.logger.add()  # write data to disc to be plotted
    es.disp()

# Display results
es.result_pretty()

# 2D Plot
fig, ax = plt.subplots()
ax.contour(X, Y, Z, 40, cmap=cm.coolwarm)
ax.axis('equal')
ax.scatter( es.mean[0], es.mean[1], marker='o', color='black', s=50, zorder=10 )

es.result_pretty()

cma.plot()

optimize and fmin2 interfaces

These are more simple interfaces that give little control over particular details of the looping, but work the same for optimization.

es = cma.CMAEvolutionStrategy(8 * [0], 0.5)
es.optimize( f )

xopt, es = cma.fmin2( f, 8 * [0], 0.5)

help( cma )