Note

Go to the end to download the full example code.

Optimization of the Rastrigin test function

The Rastrigin function is defined by:

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

where A=10 and \vect{x}\in[-5.12,5.12]^n.

It has a global minimum at \vect{x} = \vect{0} where f(\vect{x})= - 10.

This function has many local minima, so optimization algorithms must be run from multiple starting points.

In our example, we consider the bidimensional case, i.e. n=2.

References:

  • Rastrigin, L. A. “Systems of extremal control.” Mir, Moscow (1974).

  • Rudolph, G. “Globale Optimierung mit parallelen Evolutionsstrategien”. Diplomarbeit. Department of Computer Science, University of Dortmund, July 1990.

Definition of the problem

import openturns as ot
import openturns.viewer as otv
import math as m

ot.Log.Show(ot.Log.NONE)


def rastriginPy(X):
    A = 10.0
    delta = [x**2 - A * m.cos(2 * m.pi * x) for x in X]
    y = A + sum(delta)
    return [y]


dim = 2
rastrigin = ot.PythonFunction(dim, 1, rastriginPy)
print(rastrigin([1.0, 1.0]))

[-8]

Making rastrigin into a MemoizeFunction will make it recall all evaluated points.

rastrigin = ot.MemoizeFunction(rastrigin)

This example is academic and the point achieving the global minimum of the function is known.

xexact = [0.0] * dim
print(xexact)

[0.0, 0.0]

The optimization bounds must be specified.

lowerbound = [-4.4] * dim
upperbound = [5.12] * dim
bounds = ot.Interval(lowerbound, upperbound)

Plot the iso-values of the objective function

graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
graph.setLegendPosition("upper left")
graph.setLegendCorner([1.0, 1.0])
view = otv.View(graph)
Rastrigin function

We see that the Rastrigin function has several local minima. However, there is only one single global minimum at \vect{x}^\star=(0, 0).

Create the problem and set the optimization algorithm

problem = ot.OptimizationProblem(rastrigin)

We use the Cobyla algorithm and run it from multiple starting points selected by a LowDiscrepancyExperiment.

size = 64
distribution = ot.JointDistribution([ot.Uniform(lowerbound[0], upperbound[0])] * dim)
experiment = ot.LowDiscrepancyExperiment(ot.SobolSequence(), distribution, size)
solver = ot.MultiStart(ot.Cobyla(problem), experiment.generate())

Visualize the starting points of the optimization algorithm

startingPoints = solver.getStartingSample()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(startingPoints)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegends([""])
# sphinx_gallery_thumbnail_number = 2
view = otv.View(graph)
Rastrigin function

We see that the starting points are well spread across the input domain of the function.

Solve the optimization problem

solver.run()
result = solver.getResult()
xoptim = result.getOptimalPoint()
print(xoptim)

[1.60785e-06,-3.76144e-06]

xexact

[0.0, 0.0]

We can see that the solver found a very accurate approximation of the exact solution.

Analyze the optimization process

MultiStart ran an instance of Cobyla from each starting point.

Let us focus on the instance that found the global minimum. How many times did it evaluate rastrigin?

result.getCallsNumber()

2722

Let us view these evaluation points.

inputSample = result.getInputSample()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(inputSample)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegendCorner([1.0, 1.0])
graph.setLegendPosition("upper left")
view = otv.View(graph)
Rastrigin function

How fast did it find the global minimum?

graph = result.drawOptimalValueHistory()
view = otv.View(graph)
Optimal value history

Let us now analyze the MultiStart process as a whole.

Since rastrigin is a MemoizeFunction, it has a getInputHistory() method which lets us see all points it was evaluated on since its creation.

inputSample = rastrigin.getInputHistory()
graph = rastrigin.draw(lowerbound, upperbound, [100] * dim)
graph.setTitle("Rastrigin function")
cloud = ot.Cloud(inputSample)
cloud.setPointStyle("bullet")
cloud.setColor("black")
graph.add(cloud)
graph.setLegendCorner([1.0, 1.0])
graph.setLegendPosition("upper left")
view = otv.View(graph)
Rastrigin function

How many times did all Cobyla instances combined call rastrigin?

print(rastrigin.getInputHistory().getSize())

2722

otv.View.ShowAll()