Algorithmes d’optimisation 1D

Objectifs de ce TP

• Comment optimiser une fonction en dimension 1 ? Comment fonctionnent les méthodes de la dichotomie et du nombre d’or ?
• Utiliser les solveur de bases de scipy comme minimize, bisect, etc. du module optimize.

Méthode de la dichotomie

Dans cette partie, nous allons chercher une racine d’une fonction f continue sur un intervalle [a,b] \subset \mathbb{R} telle que f(a)f(b)\leq 0 en utilisant la méthode de la dichotomie¹(🇬🇧: bisection method).

¹ Voir : https://fr.wikipedia.org/wiki/M%C3%A9thode_de_dichotomie

La dichotomie est une méthode d’ordre 0, c’est à dire qui n’utilise que des évaluations de la fonction f. Dans les prochains TPs, nous introduirons des méthodes d’ordre supérieur (1 pour l’évaluation de la fonction et son gradient, 2 pour la hessienne).

    \begin{algorithm}
    \caption{Dichotomie pour trouver un zéro d'une fonction (i.e., $x$ t.q. $f(x)\approx 0$)}
    \begin{algorithmic}
    \INPUT $f, a, b, \epsilon$
    \STATE $x^{\downarrow}=\min(a,b),\quad x^{\uparrow}=\max(a,b)$
    \COMMENT{Initialisation}
    \WHILE{$x^{\uparrow}-x^{\downarrow}>\epsilon$}
        \STATE $\bar{x} \gets \frac{1}{2}(x^{\uparrow}+x^{\downarrow})$
        \COMMENT{Calcul du milieu du segment courant}
    \IF{$f(x^{\downarrow}) \cdot f(\bar{x}) \leq 0$ }
        \COMMENT{Il y a une racine entre $\bar{x}$ et $x^{\uparrow}$}
        \STATE $x^{\uparrow}\gets \bar{x}$
    \ELSE
    \STATE $x^{\downarrow}\gets \bar x$
    \COMMENT{Il y a une racine entre $x^{\downarrow}$ et $\bar{x}$}
    \ENDIF
    \ENDWHILE
    \OUTPUT $\bar x$
    \end{algorithmic}
    \end{algorithm}

Question : Premier algorithme d’optimisation

Implémenter l’algorithme de la dichotomie. Le tester sur la fonction f_{test}:x\mapsto x-1 sur l’intervalle [-2,2] (on a bien f_{test}(-2)f_{test}(2)\leq 0).

La dichotomie trouve la racine d’une fonction, est-ce que l’on peut l’utiliser pour trouver le minimum de la fonction f:x\mapsto x^2+1 ?

Question : Arrêter un algorithme

Modifier l’algorithme pour ne plus utiliser une boucle while mais une boucle for avec un nombre d’itération maximum et ainsi éviter une boucle infinie. On pourra arrêter l’algorithme avant d’avoir atteint le nombre maximum d’itérations en utilisant un break².

² Voir : https://docs.python.org/3/tutorial/controlflow.html#break-and-continue-statements-and-else-clauses-on-loops

for i in range(5):
  if i == 3:
    break
  print(i)

Le critère d’arrêt peut aussi être modifié. On peut penser à contrôler la taille de l’intervalle courant, la norme L^2 entre les bornes de l’intervalle, la distance pour l’itéré courant à 0 en norme L^p

Question : Application de la dichotomie

Utiliser cet algorithme pour trouver le minimum global de la fonction suivante x\mapsto x^4 - (x-1)^2+1. On rappelle que trouver un minimum d’une fonction revient à annuler son gradient (sa dérivée en dimension 1).

Pour visualiser la convergence, on tracera la taille de l’intervalle courant en fonction de l’indice de l’itération (on pourra modifier l’algorithme pour sortir la liste des solutions courantes, en utilisant append pour ajouter un élément à une liste. Un exemple d’utilisation de cette fonction est proposé ci-dessous). Quelle est l’influence de la précision du critère de convergence sur le nombre d’itérations effectuées ? Comparer vos résultats avec l’implémentation de scipy.optimize.bisect.

my_list = []
my_list.append(3)
print("Ma liste est:", my_list)
my_list.append([5,6])
print("Ma liste est:", my_list)
print("Taille de la liste:", len(my_list))

Ma liste est: [3]
Ma liste est: [3, [5, 6]]
Taille de la liste: 2

Pour aller plus loin

En pratique il peut être plus informatif de visualiser la convergence en fonction du temps de calcul plutôt qu’en fonction du nombre d’itération. C’est notamment le cas si l’on compare entre eux deux algorithmes. Les itérations ont des coûts de calcul très différents. Proposez une telle visualisation, par exemple avec³ un code similaire à :

    import time
    t0 = time.perf_counter()
    ... # Instruction dont on veut mesurer le temps d'exécution
    t1 = time.perf_counter()
    print(t1-t0)

5.042197881266475e-05

D’autres visualisations sont possibles, voir par exemple les graphiques générés par benchopt⁴.

⁴ Voir : https://benchopt.github.io/results/

³ Voir : https://docs.python.org/3/library/time.html

Recherche de minimum local par la méthode du nombre d’or

La méthode du nombre d’or (🇬🇧: golden-section search) permet de trouver le minimum d’une fonction f continue et unimodale sur l’intervalle [a,b]\subset\mathbb{R}. On note par la suite \varphi=\frac{1+\sqrt{5}}{2} le nombre d’or. Cette méthode est illustrée dans la figure ci-dessous et décrite par l’algorithme suivant :

Illustration des grandeurs utiles dans la méthode du nombre d’or

    \begin{algorithm}
    \caption{La méthode du nombre d'or}
    \begin{algorithmic}
    \INPUT $f, a, b, \epsilon$
    \COMMENT{Initialisation:}
    \STATE $x^{\downarrow}=\min(a,b), x^{\uparrow}=\max(a,b)$
    \STATE $x_1 \gets \frac{1}{\varphi} x^{\downarrow} + \left(1-\frac{1}{\varphi}\right)x^{\uparrow}$
    \STATE $f_1 \gets f(x_1)$
    \STATE $x_2 \gets \frac{1}{\varphi} x^{\uparrow} + \left(1-\frac{1}{\varphi}\right)x^{\downarrow}$
    \STATE $f_2 \gets f(x_2)$
    \WHILE{$x^{\uparrow}-x^{\downarrow}>\epsilon$}
    \IF{$f_1 < f_2$}
        \STATE $x^\uparrow \gets x_2$
        \STATE $x_2 \gets x_1$
        \STATE $f_2 \gets f_1$
        \STATE $x_1 \gets \frac{1}{\varphi} x^{\downarrow} + \left(1-\frac{1}{\varphi}\right)x^{\uparrow}$
        \STATE $f_1 \gets f(x_1)$
    \ELSE
    \STATE $x^\downarrow \gets x_1$
    \STATE $x_1 \gets x_2$
    \STATE $f_1 \gets f_2$
    \STATE $ x_2 \gets \frac{1}{\varphi} x^{\uparrow} + \left(1-\frac{1}{\varphi}\right)x^{\downarrow}$
    \STATE $ f_2 \gets f(x_2) $

    \ENDIF
    \ENDWHILE
    \OUTPUT $\frac{x^\uparrow + x^\downarrow}{2}$
    \end{algorithmic}
    \end{algorithm}

Question : Méthode du nombre d’or et extremum

Implémenter la méthode du nombre d’or pour une fonction f quelconque en entrée, en utilisant une boucle for et un break.

Question : Applications

Tester la méthode du nombre d’or sur différentes fonctions :

• x\mapsto x^2+2 sur [-5, 5] puis sur [-5, 20],
• x\mapsto x^6 + 3 \exp(-x^2) + \frac{1}{2}\sin\left(\frac{5x}{2}\right) sur [0.3, 1.5].

On fera une vérification visuelle ici des hypothèses. Comparer votre code avec l’implémentation de scipy.optimize.golden.

Question : Influence de l’initialisation

Reprendre la dernière fonction sur l’intervalle \left[-\frac{3}{2}, \frac{3}{2} \right] en prenant différentes initialisations plusieurs fois. Que remarquez-vous?

Pour aller plus loin

Chaque algorithme a des hypothèses différentes et tous ne vont pas forcément converger vers le même minimum suivant les initialisations. Pour vérifier que nos algorithmes convergent bien, on peut utiliser des méthodes déjà implémentées.

Optimisation avec scipy

Le package scipy a notamment une fonction minimize dans son module optimize⁵. Afficher la documentation de cette fonction dans le programme python à l’aide de la commande help (on pourra utiliser le raccourci ?).

⁵ Voir : https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

Question : Avec scipy

Regarder attentivement :

• les méthodes disponibles,
• l’initialisation (en déduire un avantage de ces méthodes par rapport à la dichotomie),
• la tolérance,
• le rappel (🇬🇧: callback).

Reprendre la deuxième fonction de la question précédente et utiliser différents algorithmes disponibles, on comparera les résultats et les temps de calcul pour une tolérance de 10^{-1}, 10^{-5} puis 10^{-9} avec la même initialisation (les temps de calculs pourront être moyennés sur plusieurs répétitions pour plus de précision). Où se situe la dichotomie par rapport aux autres algorithmes ?

Question : Pour aller plus loin avec les callback

Si vous étudiez bien la documentation de cette fonction minimize, l’utilisateur n’a pas d’accès direct au numéro de l’itération (seulement aux itérés générés). Utilisez le callback (et args si vous le souhaitez) pour afficher le numéro de l’itération k, le point x_k considéré et la valeur de f(x_k).

# Exemple d'utilisation de callback
iterations = []  # liste pour stocker les itérés

def callback(x):  # callback ne prend que l'itéré courant
    iterations.append(x)  # on ajoute l'itéré courant

scipy.optimize.minimize(ma_fonction, x0, callback=callback)
print(iterations)  # liste des itérés

Stockez ces résultats et visualisez la trajectoire engendrée par (f(x_k))_k.

Dans la dernière fonction considérée, l’hypothèse d’unimodalité n’était pas vérifiée sur l’intervalle considéré. Les méthodes d’optimisation en général ne vont pas trouver le minimum global, mais un minimum local. Sortir des minima locaux est un problème qui requiert des techniques plus avancées⁶.

⁶ Par exemple: https://fr.wikipedia.org/wiki/Algorithme_du_gradient_stochastique

Question : grille de points d’initialisation

Utilisez la méthode par défaut de scipy avec une grille de points d’initialisation sur la fonction g: x \mapsto \exp\left(\frac{1}{2}x^2\right) - x^3 + 10\cos(x), pour x\in [-3, 3]. On parle de bassin d’attraction pour un ensemble de point qui convergent vers un même minimum. Combien de bassin d’attractions observez-vous ? Associez leur une couleur et proposez une visualisation de la correspondance entre le point initial et la solution trouvée.