.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_numerical_methods_optimization_plot_optimization_rosenbrock.py:
Quick start guide to optimization
=================================
In this example, we perform the optimization of the Rosenbrock test function.
Let :math:`a, b\in\mathbb{R}` be parameters. The Rosenbrock function is defined by
.. math::
f(x_1, x_2) = (a-x_1)^2 + b(x_2 - x_1^2)^2
for any :math:`\mathbf{x}\in\mathbb{R}^2`.
This function is often used with :math:`a=1` and :math:`b=100`. In this case, the function has a single global minimum at:
.. math::
\mathbf{x}^\star = (1, 1)^T.
This function has a nonlinear least squares structure.
*References*
* Rosenbrock, H.H. (1960). "An automatic method for finding the greatest or least value of a function". The Computer Journal. 3 (3): 175–184.
Definition of the function
--------------------------
.. code-block:: default
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
.. code-block:: default
rosenbrock = ot.SymbolicFunction(['x1', 'x2'],
['(1-x1)^2+100*(x2-x1^2)^2'])
.. code-block:: default
x0 = ot.Point([-1.0, 1.0])
.. code-block:: default
xexact = ot.Point([1.0, 1.0])
.. code-block:: default
lowerbound = [-2.0, -2.0]
upperbound = [2.0, 2.0]
Plot the iso-values of the objective function
---------------------------------------------
.. code-block:: default
rosenbrock = ot.MemoizeFunction(rosenbrock)
.. code-block:: default
graph = rosenbrock.draw(lowerbound, upperbound, [100]*2)
graph.setTitle("Rosenbrock function")
view = viewer.View(graph)
.. image:: /auto_numerical_methods/optimization/images/sphx_glr_plot_optimization_rosenbrock_001.png
:alt: Rosenbrock function
:class: sphx-glr-single-img
We see that the minimum is on the top right of the picture and the starting point is on the top left of the picture. Since the function has a long valley following the curve :math:`x_2 - x^2=0`, the algorithm generally have to follow the bottom of the valley.
Create and solve the optimization problem
-----------------------------------------
.. code-block:: default
problem = ot.OptimizationProblem(rosenbrock)
.. code-block:: default
algo = ot.Cobyla(problem)
algo.setMaximumRelativeError(1.e-1) # on x
algo.setMaximumEvaluationNumber(50000)
algo.setStartingPoint(x0)
algo.run()
.. code-block:: default
result = algo.getResult()
.. code-block:: default
xoptim = result.getOptimalPoint()
xoptim
.. raw:: html
[0.99251,0.985022]
.. code-block:: default
delta = xexact - xoptim
absoluteError = delta.norm()
absoluteError
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.01674594609725928
We see that the algorithm found an accurate approximation of the solution.
.. code-block:: default
result.getOptimalValue() # f(x*)
.. raw:: html
[5.6392e-05]
.. code-block:: default
result.getEvaluationNumber()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
10520
.. code-block:: default
graph = rosenbrock.draw(lowerbound, upperbound, [100]*2)
cloud = ot.Cloud(ot.Sample([x0, xoptim]))
cloud.setColor("black")
cloud.setPointStyle("bullet")
graph.add(cloud)
graph.setTitle("Rosenbrock function")
view = viewer.View(graph)
.. image:: /auto_numerical_methods/optimization/images/sphx_glr_plot_optimization_rosenbrock_002.png
:alt: Rosenbrock function
:class: sphx-glr-single-img
We see that the algorithm had to start from the top left of the banana and go to the top right.
.. code-block:: default
graph = result.drawOptimalValueHistory()
view = viewer.View(graph)
.. image:: /auto_numerical_methods/optimization/images/sphx_glr_plot_optimization_rosenbrock_003.png
:alt: Optimal value history
:class: sphx-glr-single-img
The function value history make the path of the algorithm clear. In the first step, the algorithm went in the valley, which made the function value decrease rapidly. Once there, the algorithm had to follow the bottom of the valley so that the function decreased but slowly. In the final steps, the algorithm found the neighbourhood of the minimum so that the local convergence could take place.
In order to see where the function was evaluated, we use the `getInputSample` method.
.. code-block:: default
inputSample = result.getInputSample()
.. code-block:: default
graph = rosenbrock.draw(lowerbound, upperbound, [100]*2)
graph.setTitle("Rosenbrock function solved with Cobyla")
cloud = ot.Cloud(inputSample)
graph.add(cloud)
view = viewer.View(graph)
.. image:: /auto_numerical_methods/optimization/images/sphx_glr_plot_optimization_rosenbrock_004.png
:alt: Rosenbrock function solved with Cobyla
:class: sphx-glr-single-img
We see that the algorithm made lots of evaluations in the bottom of the valley before getting in the neighbourhood of the minimum.
Solving the problem with NLopt
------------------------------
We see that the `Cobyla` algorithm required lots of function evaluations. This is why we now use the `NLopt` class with the LBFGS algorithm. However, the algorithm may use input points which are far away from the input domain we used so far. This is why we had bounds to the problem, so that the algorithm never goes to far away from the valley.
.. code-block:: default
bounds = ot.Interval(lowerbound, upperbound)
.. code-block:: default
problem = ot.OptimizationProblem(rosenbrock)
problem.setBounds(bounds)
.. code-block:: default
algo = ot.NLopt(problem, 'LD_LBFGS')
algo.setStartingPoint(x0)
algo.run()
.. code-block:: default
result = algo.getResult()
.. code-block:: default
xoptim = result.getOptimalPoint()
xoptim
.. raw:: html
[1,1]
.. code-block:: default
delta = xexact - xoptim
absoluteError = delta.norm()
absoluteError
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
7.740583643426769e-12
We see that the algorithm found an extremely accurate approximation of the solution.
.. code-block:: default
result.getOptimalValue() # f(x*)
.. raw:: html
[1.77616e-23]
.. code-block:: default
result.getEvaluationNumber()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
44
This number of iterations is much less than the previous experiment.
.. code-block:: default
inputSample = result.getInputSample()
.. code-block:: default
graph = rosenbrock.draw(lowerbound, upperbound, [100]*2)
graph.setTitle("Rosenbrock function solved with NLopt/LD_LBFGS")
cloud = ot.Cloud(inputSample)
graph.add(cloud)
view = viewer.View(graph)
plt.show()
.. image:: /auto_numerical_methods/optimization/images/sphx_glr_plot_optimization_rosenbrock_005.png
:alt: Rosenbrock function solved with NLopt/LD_LBFGS
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.649 seconds)
.. _sphx_glr_download_auto_numerical_methods_optimization_plot_optimization_rosenbrock.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_optimization_rosenbrock.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_optimization_rosenbrock.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_