In [1]:
import numpy as np
import numpy.linalg as lin
import matplotlib.pylab as plt
import plotly.offline as py
import plotly.graph_objects as go

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

np.set_printoptions(precision=3, linewidth=150, suppress=True)
plt.style.use(['seaborn-whitegrid','data/cours.mplstyle'])

Problème d'optimisation¶

Un problème d'optimisation se présente sous la forme suivante :

Soit une fonction $J : ℝ^n \rightarrow ℝ$, trouver le minimum de $J$ c.a.d. trouver ${\bf u}$ tel que

$$J({\bf u}) = \inf_{{\bf v} \in ℝ^n} J({\bf v})$$

Problème d'optimisation avec contrainte¶

Il est possible de chercher u non pas dans $ℝ^n$ mais dans une partie de $ℝ^n$, il s'agit alors d'un problème d'optimisation avec contrainte.

Exemple On est en 2D et on cherche le minimum de $J(x,y)$ mais avec la contrainte que $y > x$. Cela revient à chercher dans la partie de $ℝ^2$ qui vérifie $y > x$.

Nous ne regarderons pas les problèmes d'optimisation avec contraintes dans ce cours mais il s'agit de problèmes importants que vous verrez dans votre scolarité.

La méthode du gradient¶

Une méthode pour résoudre un problème d'optimisation en 2D est de l'imaginer comme un terrain avec du relief. La fonction $J(x,y)$ représente l'altitude en tout point $(x,y)$.

Pour trouver le minimum de $J$ il suffit de partir d'un point au hasard et de descendre dans la direction qui descend le plus. Une goutte d'eau suit ce chemin.

C'est la méthode du gradient.

Pour bien comprendre la méthode il est nécessaire de bien comprendre les dérivées partielles ainsi que les notations qui sont résumée dans ce mémo sur les dérivées partielles.

L'algorithme du gradient est donc :

  • prendre un point de départ au hasard ${\bf p^0} = (x_0, y_0)$
  • calculer le gradient de J en ce point

$$ \nabla J(x_0, y_0) = \begin{bmatrix} \frac{\partial J}{\partial x} \\ \frac{\partial J}{\partial y} \end{bmatrix} (x_0, y_0) $$

  • avancer dans la direction opposée (le grandiant monte) : ${\bf p}^{k+1} = {\bf p}^k - \mu \, \nabla J({\bf p}^k)$

et on recommence cette dernière étape jusqu'à ce qu'on arrive à un point fixe c.a.d. que $|| {\bf p}^{k+1} - {\bf p}^k|| < \varepsilon$ avec $\varepsilon$ une toute petite valeur.

In [2]:
def J(x, y):
    return x**2 + 0.5 * y**2  - 2 * x + 3
In [3]:
x = np.linspace(-3,3,100)
y = np.linspace(-3,3,100)

mx, my = np.meshgrid(x,y)
mz = J(mx, my)

trace = go.Surface(x=mx, y=my, z=mz)
layout = go.Layout(scene = {'aspectmode':'data'})
fig = go.Figure(data=[trace], layout=layout)
py.iplot(fig, show_link=False)
In [4]:
def grad_J(x,y):
    return np.array([2*x-2, y])   # calculé à la main à partir de J (on aimerait que cela soit automatique)

Si vous ne comprennez pas le calcul du gradiant de $J$, n'allez pas plus loin. Retrouvez votre cours sur ce sujet, regardez le mémo, posez des questions.

In [5]:
# Algo du gradient pour trouver le minimum

x = np.array([0,0])  # un point au hasard, changez pour voir
µ = 0.1     # plus il est petit et moins on avance vite, cf formule ci-dessus
e = 0.0001  # mon epsilon pour la condition d'arrêt

while True:
    x_old = x
    x =  x - µ * grad_J(*x)  # *x donne en arguments toutes les valeurs de x donc x[0] en 1er arg et x[1] en 2e
    if np.square(x_old - x).sum() < e**2:
        break
print("Le minimum est au point ", x, "La valeur de J en ce point est ", J(*x))

Le minimum est au point  [1. 0.] La valeur de J en ce point est  2.0000001053122918

Je vous invite à bouger la figure pour être au dessus, zoomer pour voir ce qui va se passer et maintenant déplacer la souris sur la figure. Vous devez voir une ellipse qui montre la ligne de niveau (tous les points ayant la même valeur de z). Le minimum est lorsque cette ellipse n'est qu'un point. C'est justement en [1, 0, 2].

Étude de la convergence du gradient¶

On stocke toutes les valeurs des points entre notre point initial et la solution pour pouvoir tracer des courbes de convergence.

Donc le résultat de notre méthode de gradient doit être un ensemble de points.

In [6]:
def minimum_J(start_value, µ=0.1, e = 0.001):
    x = [np.array(start_value)]
    while True:
        x.append(x[-1] - µ * grad_J(*x[-1]))
        if np.square(x[-1] - x[-2]).sum() < e**2:
            break
        # la suite n'est que des tests pour se protéger
        if np.square(x[-1] - x[-2]).sum() > 1E6:  # au cas où on diverge
            print("DIVERGE")
            break
        if len(x) > 200:  # c'est trop long, je crains la boucle infinie
            print('Trop long, boucle infinie ?')
            break
    return np.array(x)


x = minimum_J(start_value = (0,1))  # je prends une valeur initiale qui n'est pas alignée avec la solution
In [7]:
plt.plot(x[:,0], x[:,1], 'x:')
Out[7]:
[<matplotlib.lines.Line2D at 0x7fbc5e968c50>]
No description has been provided for this image
In [8]:
fig = go.Figure(data=[go.Scatter3d(x=x[:,0], y=x[:,1], z=J(x[:,0], x[:,1]), marker={'size':3})])
fig.show()

C'est joli, n'est-ce pas ?

Impact de µ¶

Regardons comment µ influence sur la convergence c.a.d. une fois qu'on a la direction de la plus grande pente en un point de combien doit-on avancer. Avec un µ grand on fait des grands pas, avec un µ petit on fait des petits pas.

Avec µ = 0.1 on a fait des petits pas.

In [9]:
# µ = 2

x = minimum_J(start_value = (0,1), µ = 2)
plt.plot(x[:,0], x[:,1], 'x:')

DIVERGE
Out[9]:
[<matplotlib.lines.Line2D at 0x7fbc5c93f9b0>]
No description has been provided for this image
In [10]:
fig = go.Figure(data=[go.Scatter3d(x=x[:,0], y=x[:,1], z=J(x[:,0], x[:,1]), marker={'size':3})])
fig.show()

Pas bon, on a divergé. Les pas sont trop grands. On va dans la bonne direction mais tellement loin qu'on remonte de l'autre coté de la cuvette.

On voit bien qu'en y on passe de 1 a -1 puis à 1, -1 etc sans jamais s'arrêter à 0 qui est la solution pour y.

Question Pourquoi la valeur de x explose au lieu de faire comme y ?

In [11]:
# µ = 0.8

x = minimum_J(start_value = (0,1), µ = 0.8)
plt.plot(x[:,0], x[:,1], 'x:')
Out[11]:
[<matplotlib.lines.Line2D at 0x7fbc5c54e080>]
No description has been provided for this image
In [12]:
fig = go.Figure(data=[go.Scatter3d(x=x[:,0], y=x[:,1], z=J(x[:,0], x[:,1]), marker={'size':3})])
fig.show()

On a convergé mais c'est moins joli que la première fois. Par contre

In [13]:
len(x)
Out[13]:
17

on a convergé en 17 itérations. Regardons combien il a fallu d'itération pour µ = 0.1 :

In [14]:
x = minimum_J(start_value = (0,1), µ = 0.1)
len(x)
Out[14]:
46

µ = 0.8 donne une courbe moins jolie mais on converge presque 3 fois plus vite qu'avec µ = 0.1. Donc c'est mieux.

In [15]:
# µ = 1

x = minimum_J(start_value = (0,1), µ = 1)
plt.plot(x[:,0], x[:,1], 'x:')

Trop long, boucle infinie ?
Out[15]:
[<matplotlib.lines.Line2D at 0x7fbc5ec1cc88>]
No description has been provided for this image
In [16]:
fig = go.Figure(data=[go.Scatter3d(x=x[:,0], y=x[:,1], z=J(x[:,0], x[:,1]), marker={'size':3})])
fig.show()

On est parti de [0,1] pour aller en [2,0] puis en [0,0] puis en [2.0].

Comme on a vu que depuis [2,0] on va en [0,0] cela veut dire qu'on oscille infiniment entre [0,0] et [2,0].

Sachant que la solution recherchée est [1,0] cela veut dire qu'on va bien dans la bonne direction mais on fait un pas 2 fois trop grand à chaque fois.

Moralité La valeur de µ est importante. Si elle est trop petite on perd du temps, si elle est trop grande on ne trouve pas la solution.

In [ ]: