Pour ce TP noté, vous devrez charger votre fichier sur Moodle, avant 18h, mercredi 26/04/2023.
Important
Pour ce travail vous devez déposer un unique fichier au format .py, structuré avec des cellules, dont le nom est examen_tp_HAX606X_Prenom_Nom.py. Vous remplirez votre nom, prénom de manière adéquate. Si le nom de fichier ne suit pas cette forme, vous perdez 1 pt.
De plus, les personnes qui n’auront pas soumis leur devoir sur Moodle avant la limite obtiendront comme note zéro. Ainsi, pensez à soumettre une version provisoire par sécurité à la fin de la séance ! (aucune exception ne sera faite)
La note totale est sur 20 points répartis comme suit :
qualité des réponses aux questions : 15 pts,
qualité de rédaction et d’orthographe : 1 pts,
qualité des graphiques (légendes, couleurs) : 1 pt
qualité du code (noms de variables clairs, commentaires, code synthétique, etc.) : 2 pt
reproductibilité (exécution correcte sur la machine du correcteur) et absence de bug : 1 pt
La descente de gradient classique permet de minimiser une fonction convexe f: \mathbb{R}^d \to \mathbb{R}:
x^\star \in \argmin_{x\in \mathbb{R}^d} f(x) \enspace.
Remarque: on note \in \argmin pour désigner une solution du problème d’optimisation, le minimiseur n’étant pas forcément unique.
Pour ce faire, l’étape de descente s’effectue à l’aide de l’équation de mise à jour
x_{k+1} \gets x_{k} - \gamma \nabla f(x_k) \enspace,
où \gamma>0 est le pas de descente.
Il est cependant fréquent de vouloir minimiser une fonction avec des contraintes. Ici, on représente les contraintes par un ensemble \mathcal{C} \subseteq \mathbb{R}^d, où \mathcal{C} est un ensemble fermé et convexe. Le problème d’optimisation contraint s’écrit alors:
x^\star \in \argmin_{x\in \mathcal{C}} f(x) \enspace.
\tag{1} Une stratégie naturelle est d’utiliser la descente de gradient projeté pour résoudre ce problème:
x_{k+1} \gets \Pi_\mathcal{\mathcal{C}}\left(x_{k} - \gamma \nabla f(x_k)\right) \enspace,
où la fonction \Pi_\mathcal{\mathcal{C}} est la projection sur l’ensemble \mathcal{C}.
Question 1 : Visualisation des contraintes et des lignes de niveaux
Cherchons à minimiser la fonction quadratique suivante sur \mathcal{C}=[0, 1]^2 (ici d=2) :
f_{test}:x\mapsto \frac{1}{2}(x-x_0)^\top Q (x-x_0)\enspace,
avec x_0=(1.5, 0)^\top et Q=\begin{pmatrix} 3 & 2 \\ 2 & 3\end{pmatrix} une matrice définie positive.
Écrire un code permettant de reproduire la Figure 1 ci-dessous. On pourra utiliser l’outils mpatches venant de matplotlib.patches pour délimiter la zone de contrainte (ou proposer une autre visualisation).
Figure 1: Fonction quadratique à minimiser et contraintes
Question 2 : Descente de gradient et descente de gradient projeté
Coder la descente de gradient classique (🇬🇧: gradient descent (GD)) avec un critère d’arrêt sur la norme du gradient courant. Prendre pour pas de descente \gamma=0.05 et tracer la trajectoire des itérés avec en point initial x_{init}=(0.2, 0.5)^\top. Faire de même pour la descente de gradient projeté et visualiser pour chaque méthode la trajectoire des itérés obtenue (on pourra s’inspirer de l’ex)
Figure 2: Trajectoires des itérés pour la descente de gradient et la descente de gradient projeté
Une autre méthode pour l’optimisation sous contrainte: la descente miroir
La descente de gradient projeté (🇬🇧: projected gradient descent, PGD) suit les itérés de la descente de gradient classique puis brutalement contraint sa trajectoire à rester dans la zone de contraintes. Essayons de mitiger ce comportement brutal de la descente de gradient projeté avec la descente miroir (🇬🇧: mirror descent, MD).
Tout d’abord, notons qu’une étape de descente de gradient projeté peut s’écrire:
x_{k+1} \gets \argmin_{x\in\mathcal{C}} \langle \nabla f(x_k),x\rangle + \frac{1}{2\gamma}\| x-x_k\|^2_2\enspace.
\tag{2} Cette formulation fait intervenir la norme \ell_2 (ou norme euclidienne standard) sur \mathbb{R}^d. Mais on peut généraliser cette formulation en prenant non plus la norme \ell_2 mais une divergence de Bregman:
D_\Psi(x, y) = \Psi(x) - \Psi(y) - \langle \nabla \Psi(y), x-y\rangle\enspace,
\tag{3} pour un \Psi: \mathcal{C} \to \mathbb{R} différentiable sur l’intérieur de \mathcal{C} .
La descente miroir est alors définie par:
x_{k+1} \gets \argmin_{x\in\mathcal{C}} \langle \nabla f(x_k),x\rangle + \frac{1}{\gamma}D_\Psi(x, x_k)\enspace.
\tag{4}
Cas particulier
En choisissant \Psi(x)=\frac{1}{2}\|x\|^2_2 dans l’équation Équation 4, on retrouve la descente de gradient projeté classique, telle que décrite dans l’équation Équation 2.
Dans la suite, on suppose de plus que la fonction dite “miroir” \phi = \nabla \Psi:\mathcal{C} \to \mathbb{R}^d est bijective de \mathcal{C} vers \phi(\mathcal{C})=\mathbb{R}^d, et on admet que l’on peut écrire la descente miroir sous la forme:
x_{k+1}\gets \phi^{-1}\left(\phi(x_k) - \gamma \nabla f(x_k)\right) \enspace.
\tag{5} On visualise l’algorithme avec la Figure 3 ci-dessous:
Figure 3: Schéma de la descente miroir
Question 3 : Entropie et divergence de Bregman (d=1)
Dans le cas d=1, et \mathcal{C}=[0,1] choisissons l’entropie binaire négative définie sur ]0,1[ par \Psi(x) = x\log(x) + (1-x)\log(1-x)\enspace.
Tracer la fonction \Psi et mettre en valeur son minimum; sur des sous-graphiques juxtaposés, tracer la fonction miroir et son inverse (\phi et \phi^{-1}); on pourra s’inspirer Figure 4.
Figure 4: Visualisation de l’entropie et autres propriétés
Question 4 : Implémentation de la descente miroir (d=2)
Implémenter la descente miroir (Équation 5) avec le choix
\Psi(x)=\Psi(x_1,x_2)=\sum_{i=1}^2 x_i\log(x_i) + (1-x_i)\log(1-x_i)\enspace.
On utilisera un critère d’arrêt basé uniquement sur le nombre d’itérations pour le moment.
Tester la descente miroir sur la fonction quadratique f_{test} introduite précédemment (avec \gamma=1 par exemple) en comparant la trajectoire obtenue avec la descente de gradient projeté. Commenter le résultat obtenu.
Figure 5: Visualisation des itérés pour les trois méthodes (GD, PGD, MD).
Question 5 : Monitorer la convergence (d=2)
a) Saut de dualité
Dans cette question on va introduire le saut de dualité. C’est une fonction qui indique le niveau proximité d’un point pour le problème d’optimisation Équation 1. Ainsi, le saut de dualité tend vers zéro quand on l’évalue sur une suite de points qui converge vers la solution du problème d’optimisation Équation 1, .
Pour le définir, on introduit la fonction g:\mathbb{R^2} \to \mathbb{R} définie par
\begin{align*}
g: \mathbb{R^2} &\to \mathbb{R}\\
x & \mapsto \frac{1}{2} x^\top Q^{-1} x + x^\top x_0 \enspace,
\end{align*}
et
\begin{align*}
h: \mathbb{R^2} &\to \mathbb{R}\\
u & \mapsto h(u)=\sum_{i=1}^2 \max(0, u_i)\enspace,
\end{align*}
Enfin, en notant u(x)=Q(x-x_0), la fonction Gap est définie ainsi: \begin{align*}
Gap: [0,1]^2 &\to \mathbb{R} \\
x & \mapsto f(x) + g(u(x)) + h(-u(x)) \enspace,
\end{align*}
et pour x\notin[0,1]^2, Gap(x)=+\infty.
Afficher Gap(x_n) pour les itérés (x_{n})_{n=1,...,50} des algorithmes de descente de gradient projeté (PGD) et de descente miroir (MD). On pourra pour cela une comparaison visuelle de la convergence des deux algorithmes en affichant l’évolution de Gap(x_n) en fonction de n et du choix du pas. Ici \gamma\in\{0.1, 0.2, 0.5, 0.1, 2, 5\}. Commenter les résultats obtenus.
Figure 6: Évolution du saut de dualité avec les itérés.
b) Stabilité des itérés
A partir de la liste des itérés, on peut monitorer la stabilité de la trajectoire des itérés. Un critère d’arrêt souvent utilisé en pratique est de s’arrêter à l’étape k>1 si:
\|x_k - x_{k-1}\|_2\leq \epsilon \enspace.
Pour \gamma\in\{0.1, 0.5, 1\}, visualiser l’évolution de ce critère pour PGD et MD. Commenter le résulat obtenu (on pourra s’aider en regardant le contenu de la liste itérés obtenus).
Figure 7: Évolution de la stabilités des itérés.
Question 6 : Solution avec scipy?
Utiliser la fonction minimize de scipy pour minimiser f_{test} sur [0,1]^2. Comparer et commenter avec les résultats obtenus précédement.
Question 7 : Résolution théorique du problème
Prouvez mathétatiquement quel est le minimum global de f_{test} sur l’ensemble de contraintes \mathcal{C}=[0,1]^2.
Question 8 : Changeons l’espace de contraintes
Trouver avec scipy la solution obtenue pour le problème suivant
\argmin_{\mathcal{X}=(x_1,x_2)\in\mathbb{R}^2} f_{test}(x) \text{ sc } \begin{cases} y-x &< 0 \\ x^2 + 3 y^2 &\geq 0 \end{cases} \enspace.
Optimisation non-convexe
Dans cette partie, nous utiliserons une variation de la fonction f définie par
f(x, y) = 0.9\sin(x) \times 0.9\cos(y) + 0.1 (x^2 + y^2) - 0.5
Cette fonction n’est pas convexe, mais admet un minimum global !
Question 9 : Visualisation de la fonction
Visualiser la fonction en 3D ainsi que ses lignes de niveau.
Figure 8: Fonction et lignes de niveau.
Question 10 : Minima et points d’initialisation
Prendre 5 points uniformément au hasard avec un générateur du module np.random fixé avec une seed fixée à 1 dans le carré [-3, 3]^2 et éxécuter une descente de gradient sans contrainte avec un pas de \gamma=0.1, comparer les résultats obtenus. Commenter les minima obtenus et les problèmes rencontrés. Que faire en pratique quand on rencontre une fonction non convexe à minimiser (proposer deux solutions possibles en les commentant brièvement).
Figure 9: Initialization et convergence.
Code source
---title: "TP noté"author: "Joseph Salmon, Tanguy Lefort"toc: truelang: frformat: html: html-math-method: katex toc-location: left toc-depth: 6 code-tools: true self-contained: true page-layout: article fig-format: retina echo: false html: code-fold: true code-tools: trueexecute: warning: false error: truejupyter: python3---<!-- sources:- https://vene.ro/blog/mirror-descent.html blog Vlad Nicolae- book by Beck 2017.- https://arxiv.org/pdf/1405.4980.pdf poly Bubeck- https://tlienart.github.io/posts/2018/10/27-mirror-descent-algorithm/ blog T. Lienart- https://theses.hal.science/tel-01446492 Thèse Kwoon-->Pour ce TP noté, vous devrez charger votre fichier sur Moodle, avant **18h, mercredi 26/04/2023**.::: {.callout-important}Pour ce travail vous devez déposer un **unique** fichier au format `.py`, structuré avec des cellules, dont le nom est `examen_tp_HAX606X_Prenom_Nom.py`.Vous remplirez votre nom, prénom de manière adéquate.Si le nom de fichier ne suit pas cette forme, <span style="color:red">vous perdez 1 pt</span>.De plus, les personnes qui n'auront pas soumis leur devoir sur Moodle avant la limite obtiendront comme note <span style="color:red"> **zéro**</span>.Ainsi, pensez à soumettre une version provisoire par sécurité à la fin de la séance ! (aucune exception ne sera faite):::La note totale est sur **20** points répartis comme suit :- qualité des réponses aux questions : **15** pts,- qualité de rédaction et d’orthographe : **1** pts,- qualité des graphiques (légendes, couleurs) : **1** pt- qualité du code (noms de variables clairs, commentaires, code synthétique, etc.) : **2** pt- reproductibilité (exécution correcte sur la machine du correcteur) et absence de bug : **1** pt::: {.hidden}$$\DeclareMathOperator*{\argmin}{argmin}\DeclareMathOperator*{\argmax}{argmax}$$:::# Minimiser une fonction sous contraintes## Méthode de la descente de gradient projetéLa descente de gradient classique permet de minimiser une fonction convexe$f: \mathbb{R}^d \to \mathbb{R}$:$$x^\star \in \argmin_{x\in \mathbb{R}^d} f(x) \enspace.$$::: {.callout-note appearance="simple"}Remarque: on note $\in \argmin$ pour désigner une solution du problème d'optimisation, le minimiseur n'étant pas forcément unique.:::Pour ce faire, l'étape de descente s'effectue à l'aide de l'équation de mise à jour$$x_{k+1} \gets x_{k} - \gamma \nabla f(x_k) \enspace,$$où $\gamma>0$ est le pas de descente.Il est cependant fréquent de vouloir minimiser une fonction avec des contraintes. Ici, on représente les contraintes par un ensemble $\mathcal{C} \subseteq \mathbb{R}^d$, où $\mathcal{C}$ est un ensemble fermé et convexe.Le problème d'optimisation contraint s'écrit alors:$$x^\star \in \argmin_{x\in \mathcal{C}} f(x) \enspace.$${#eq-opt-contrainte}Une stratégie naturelle est d'utiliser la descente de gradient projeté pour résoudre ce problème:$$x_{k+1} \gets \Pi_\mathcal{\mathcal{C}}\left(x_{k} - \gamma \nabla f(x_k)\right) \enspace,$$où la fonction $\Pi_\mathcal{\mathcal{C}}$ est la projection sur l'ensemble $\mathcal{C}$.### Question 1 : Visualisation des contraintes et des lignes de niveauxCherchons à minimiser la fonction quadratique suivante sur $\mathcal{C}=[0, 1]^2$ (ici $d=2$) :$$f_{test}:x\mapsto \frac{1}{2}(x-x_0)^\top Q (x-x_0)\enspace,$$avec $x_0=(1.5, 0)^\top$ et $Q=\begin{pmatrix} 3 & 2 \\ 2 & 3\end{pmatrix}$ une matrice définie positive.Écrire un code permettant de reproduire la @fig-cons-quad ci-dessous.On pourra utiliser l'outils `mpatches` venant de `matplotlib.patches` pour délimiter la zone de contrainte (ou proposer une autre visualisation).```{python}#| code-fold: true#| fig-align: center#| label: fig-cons-quad#| fig-cap: Fonction quadratique à minimiser et contraintesimport numpy as npimport matplotlib.pyplot as pltimport matplotlib.patches as mpatchesx0 = np.array([1.5, 0.])Q = np.array([[3.0, 2.0], [2.0, 3.0]])# Define the functiondef f(x, Q=Q, x0=x0):return1/2* np.dot((x-x0).T, np.dot(Q, (x - x0)))def f_conj(x, Q=Q, x0=x0):return1/2* np.dot(x.T, np.dot(np.linalg.inv(Q), x)) + np.dot(x0.T, x)def stopping_criterion(x, x_old, eps):return np.linalg.norm(x-x_old) <= epsdef stopping_criterion_gap(x, Q, x0, eps): u = Q @ (x -x0)if np.max(x)<=1and np.min(x)>=0:return f(x, Q, x0) + f_conj(u, Q, x0) + np.sum(np.maximum(0, -u))else:return np.infx = np.linspace(-0.45, 1.55, 500)X, Y = np.meshgrid(x, x)Z = np.zeros_like(X)for i inrange(len(x)):for j inrange(len(x)): Z[i,j] = f(np.array([X[i,j], Y[i,j]]))# Create the contour plotfig, axs = plt.subplots(1,1, figsize=(5, 4))plt.contourf(X, Y, Z, cmap='RdBu_r', levels=30)axs.axes.set_aspect('equal')axs.set_xlim([-0.45, 1.55])axs.set_ylim([-0.45, 1.55])rect=mpatches.Rectangle((0,0),1,1, fill =True, color ='#FFA200', alpha=.3, linewidth =2, label="Contraintes")plt.gca().add_patch(rect)plt.title(r"$f_{test}(x)=(x-x_0)^\top Q (x-x_0)/2$")plt.legend()plt.colorbar()plt.show()```### Question 2 : Descente de gradient et descente de gradient projetéCoder la descente de gradient classique (🇬🇧: *gradient descent (GD)*) avec un critère d'arrêt sur la norme du gradient courant.Prendre pour pas de descente $\gamma=0.05$ et tracer la trajectoire des itérés avec en point initial $x_{init}=(0.2, 0.5)^\top$.Faire de même pour la descente de gradient projeté et visualiser pour chaque méthode la trajectoire des itérés obtenue (on pourra s'inspirer de l'ex)```{python}#| code-fold: true#| fig-align: center#| label: fig-trajectoires#| fig-cap: Trajectoires des itérés pour la descente de gradient et la descente de gradient projetédef projected_gd(gradf, pi, xinit, gamma, eps, maxiter): x = xinit ll = [xinit.copy()] gaps = []for i inrange(maxiter): g = gradf(x) gap = stopping_criterion_gap(x, Q, x0, eps) gaps.append(gap)if gap <= eps:breakelse: x = pi(x - gamma * g) ll.append(x)return np.array(ll), x, gapsdef proj_id(x):return xdef proj_cons(x):return np.clip(x, 0, 1)def grad_quad(x):return np.dot(Q, x-x0)ll_id, x_id, gaps_gd = projected_gd(grad_quad, proj_id, np.array([0.2, 0.5]), 0.05, 1e-7, 50)ll_pgd, x_pgd, gaps_pgd = projected_gd(grad_quad, proj_cons, np.array([0.2, 0.5]), 0.05, 1e-7, 50)init = np.array([0.2, 0.5])# x = np.linspace(-.1, 1.52, 500)X, Y = np.meshgrid(x, x)Z = np.zeros_like(X)for i inrange(len(x)):for j inrange(len(x)): Z[i,j] = f(np.array([X[i,j], Y[i,j]]))fig, axs = plt.subplots(1,2, figsize=(8, 4.6))axs[0].axes.set_aspect("equal")axs[0].contourf(X, Y, Z, cmap='RdBu_r', levels=30)rect=mpatches.Rectangle((0,0),1,1, fill=True, color='#FFA200', alpha=.3, linewidth =2, label="Contraintes")axs[0].plot(ll_id[:, 0], ll_id[:, 1], "o-", ms=5, color="red", label="GD")axs[0].plot(ll_pgd[:, 0], ll_pgd[:, 1], "o-",ms=4, color="black", label="PGD")axs[0].plot( init[0], init[1], "wo", ms=8, markeredgecolor="k", label="Initialisation")axs[0].add_patch(rect)axs[0].set_title(r"$f_{test}(x)=(x-x_0)^\top Q (x-x_0)/2$")axs[0].set_xlim([-0.45, 1.55])axs[0].set_ylim([-0.45, 1.55])# axs[0].legend()x = np.linspace(-.1, 1.52, 500)X, Y = np.meshgrid(x, x)Z = np.zeros_like(X)for i inrange(len(x)):for j inrange(len(x)): Z[i,j] = f(np.array([X[i,j], Y[i,j]]))axs[1].axes.set_aspect('equal')axs[1].contourf(X, Y, Z, cmap='RdBu_r', levels=30)rect=mpatches.Rectangle((0,0),1,1, fill=True, color='#FFA200', alpha=.3, linewidth=2, label="Contraintes")axs[1].plot(ll_id[:, 0], ll_id[:, 1], "o-", ms=5, color="red", label="GD")axs[1].plot(ll_pgd[:, 0], ll_pgd[:, 1], "o-",ms=4, color="black", label="PGD")axs[1].plot( init[0], init[1], "wo", ms=8, markeredgecolor="k", label="Initialisation")axs[1].add_patch(rect)axs[1].set_xlim([0.8, 1.2])axs[1].set_ylim([0.2, 0.6])axs[1].set_title(r"Zoom")axs[1].legend(fontsize=8)plt.show()```## Une autre méthode pour l'optimisation sous contrainte: la descente miroirLa descente de gradient projeté (🇬🇧: *projected gradient descent, PGD*) suit les itérés de la descente de gradient classique puis brutalement contraint sa trajectoire à rester dans la zone de contraintes.Essayons de mitiger ce comportement brutal de la descente de gradient projeté avec la descente miroir (🇬🇧: *mirror descent*, MD).Tout d'abord, notons qu'une étape de descente de gradient projeté peut s'écrire:$$x_{k+1} \gets \argmin_{x\in\mathcal{C}} \langle \nabla f(x_k),x\rangle + \frac{1}{2\gamma}\| x-x_k\|^2_2\enspace.$${#eq-proj-gradient}Cette formulation fait intervenir la norme $\ell_2$ (ou norme euclidienne standard) sur $\mathbb{R}^d$.Mais on peut généraliser cette formulation en prenant non plus la norme $\ell_2$ mais une [divergence de Bregman](https://en.wikipedia.org/wiki/Bregman_divergence):$$D_\Psi(x, y) = \Psi(x) - \Psi(y) - \langle \nabla \Psi(y), x-y\rangle\enspace,$${#eq-bregman}pour un $\Psi: \mathcal{C} \to \mathbb{R}$ différentiable sur l'intérieur de $\mathcal{C}$ .La descente miroir est alors définie par:$$x_{k+1} \gets \argmin_{x\in\mathcal{C}} \langle \nabla f(x_k),x\rangle + \frac{1}{\gamma}D_\Psi(x, x_k)\enspace.$${#eq-proj-bregman}::: {.callout-note}## Cas particulierEn choisissant $\Psi(x)=\frac{1}{2}\|x\|^2_2$ dans l'équation @eq-proj-bregman, on retrouve la descente de gradient projeté classique, telle que décrite dans l'équation @eq-proj-gradient.:::Dans la suite, on suppose de plus que la fonction dite "miroir" $\phi = \nabla \Psi:\mathcal{C} \to \mathbb{R}^d$ est bijective de $\mathcal{C}$ vers $\phi(\mathcal{C})=\mathbb{R}^d$, et on admet que l'on peut écrire la descente miroir sous la forme:$$x_{k+1}\gets \phi^{-1}\left(\phi(x_k) - \gamma \nabla f(x_k)\right) \enspace.$${#eq-mirror-descent}On visualise l'algorithme avec la @fig-mirror-descent ci-dessous:<!-- ::: {#fig-mirror-descent} -->{#fig-mirror-descent fig-align="center"}<!-- ::: -->### Question 3 : Entropie et divergence de Bregman ($d=1$)Dans le cas $d=1$, et $\mathcal{C}=[0,1]$ choisissons l'[entropie binaire négative](https://en.wikipedia.org/wiki/Binary_entropy_function) définie sur $]0,1[$ par$$\Psi(x) = x\log(x) + (1-x)\log(1-x)\enspace.$$Tracer la fonction $\Psi$ et mettre en valeur son minimum;sur des sous-graphiques juxtaposés, tracer la fonction miroir et son inverse ($\phi$ et $\phi^{-1}$); on pourra s'inspirer @fig-entropie.<!-- En choisissant $\Psi(x)=\frac{1}{2}\|x\|^2_2$ dans l'équation @eq-proj-bregman, on retrouve la descente de gradient projeté classique, telle que décrite dans l'équation @eq-proj-gradient. -->::: {.callout-note}## NoteLa fonction $\phi^{-1}$ obtenue ci-dessus est appelée fonction logistique et est [très utilisée en statistique et en apprentissage automatique](https://fr.wikipedia.org/wiki/R%C3%A9gression_logistique).:::```{python}#| code-fold: true#| fig-align: center#| label: fig-entropie#| fig-cap: Visualisation de l'entropie et autres propriétésdef binentropy(x):return x*np.log(x) + (1-x)*np.log(1-x)def gradbinentropy(x):return np.log(x) - np.log(1-x)def sigmoid(x):return1/(1+ np.exp(-x))x = np.linspace(0, 1, 50)y = binentropy(x)fig, axs = plt.subplots(2, 2, figsize=(6, 6))axs[0,0].plot(x, y, color="k")axs[0,0].set_ylabel(r"$\Psi(x)$")axs[0,0].set_xlabel(r"$x$")axs[0,0].plot(0.5, binentropy(0.5), "wo", color="red", ms=8, markeredgecolor="k", label="Minimum")axs[0,0].legend()axs[0,0].set_title("Entropie binaire négative")x = np.linspace(-10, 10, 101)y = sigmoid(x)axs[0,1].plot(x, y, color="k")axs[0,1].set_ylabel(r"$\phi^{-1}(y)$")axs[0,1].set_xlabel(r"$y$")axs[0,1].set_title("Fonction logistique")axs[0,1].set_xlim(-11, 11)axs[1,0].plot(y,x, color="k")axs[1,0].set_ylabel(r"$\phi(x)$")axs[1,0].set_xlabel(r"$x$")axs[1,0].set_title("Fonction miroir")axs[1,0].set_ylim(-11, 11)axs[1,1].remove()plt.tight_layout()plt.show()```### Question 4 : Implémentation de la descente miroir ($d=2$)Implémenter la descente miroir (@eq-mirror-descent) avec le choix$$\Psi(x)=\Psi(x_1,x_2)=\sum_{i=1}^2 x_i\log(x_i) + (1-x_i)\log(1-x_i)\enspace.$$On utilisera un critère d'arrêt basé uniquement sur le nombre d'itérations pour le moment.Tester la descente miroir sur la fonction quadratique $f_{test}$ introduite précédemment (avec $\gamma=1$ par exemple) en comparant la trajectoire obtenue avec la descente de gradient projeté.Commenter le résultat obtenu.```{python}#| code-fold: true#| fig-align: center#| label: fig-MD-quad#| fig-cap: Visualisation des itérés pour les trois méthodes (GD, PGD, MD).def mirror_descent(gradf, xinit, gamma, eps, maxiter): x = xinit gaps = [] ll = [xinit.copy()]for i inrange(maxiter): g = gradf(x) x = gradbinentropy(x) - gamma * g x = sigmoid(x) gap = stopping_criterion_gap(x, Q, x0, eps) gaps.append(gap)if gap <= eps:break ll.append(x)return np.array(ll), x, gapsx = np.linspace(-.1, 1.52, 500)X, Y = np.meshgrid(x, x)Z = np.zeros_like(X)for i inrange(len(x)):for j inrange(len(x)): Z[i,j] = f(np.array([X[i,j], Y[i,j]]))ll_mirror, x_mirror, gaps_mirror = mirror_descent(grad_quad, np.array([0.2, 0.5]), 1, 1e-7, 50)fig, axs = plt.subplots(1,1, figsize=(6, 4))plt.contourf(X, Y, Z, cmap='RdBu_r', levels=30)axs.axes.set_aspect('equal')rect=mpatches.Rectangle((0,0),1,1, fill=True, color='#FFA200', alpha=.3, linewidth =2, label="Contraintes")plt.plot(ll_id[:, 0], ll_id[:, 1], "o-", ms=5, color="red", label="GD")plt.plot(ll_pgd[:, 0], ll_pgd[:, 1], "o-",ms=4, color="black", label="PGD")plt.plot(ll_mirror[:, 0], ll_mirror[:, 1], "o-",ms=4, color="purple", label="MD")plt.gca().add_patch(rect)axs.plot( init[0], init[1], "wo", ms=8, markeredgecolor="k", label="Initialisation")plt.xlim([0.15, 1.15])plt.ylim([0.2, 1.2])plt.title(r"Zoom on $f_{test}(x)= (x-x_0)^\top Q (x-x_0)/2$")plt.legend()plt.show()```### Question 5 : Monitorer la convergence ($d=2$)#### a) Saut de dualitéDans cette question on va introduire le saut de dualité. C'est une fonction qui indique le niveau proximité d'un point pour le problème d'optimisation @eq-opt-contrainte. Ainsi, le saut de dualité tend vers zéro quand on l'évalue sur une suite de points qui converge vers la solution du problème d'optimisation @eq-opt-contrainte, .Pour le définir, on introduit la fonction $g:\mathbb{R^2} \to \mathbb{R}$ définie par$$\begin{align*}g: \mathbb{R^2} &\to \mathbb{R}\\ x & \mapsto \frac{1}{2} x^\top Q^{-1} x + x^\top x_0 \enspace,\end{align*}$$et$$\begin{align*}h: \mathbb{R^2} &\to \mathbb{R}\\ u & \mapsto h(u)=\sum_{i=1}^2 \max(0, u_i)\enspace,\end{align*}$$Enfin, en notant $u(x)=Q(x-x_0)$, la fonction $Gap$ est définie ainsi:$$\begin{align*}Gap: [0,1]^2 &\to \mathbb{R} \\ x & \mapsto f(x) + g(u(x)) + h(-u(x)) \enspace,\end{align*}$$et pour $x\notin[0,1]^2, Gap(x)=+\infty$.Afficher $Gap(x_n)$ pour les itérés $(x_{n})_{n=1,...,50}$ des algorithmes de descente de gradient projeté (PGD) et de descente miroir (MD).On pourra pour cela une comparaison visuelle de la convergence des deux algorithmes en affichant l'évolution de $Gap(x_n)$ en fonction de $n$ et du choix du pas. Ici $\gamma\in\{0.1, 0.2, 0.5, 0.1, 2, 5\}$.Commenter les résultats obtenus.```{python}#| code-fold: true#| fig-align: center#| label: fig-gap#| fig-cap: Évolution du saut de dualité avec les itérés.steps = np.array([0.1, 0.2, 0.5, 1, 2, 5])markers = ["o", "v", "s", "D", "X", "P"]steps_norm = np.linspace(0.1,0.95, len(steps))fig_gaps, axs_gaps = plt.subplots(1, 2, figsize=(8, 3), sharey=True)for i, step inenumerate(steps): ll_mirror, x_mirror, gaps_mirror = mirror_descent(grad_quad, np.array([0.2, 0.5]), step, 1e-7, 50) ll_pgd, x_pgd, gaps_pgd = projected_gd(grad_quad, proj_cons, np.array([0.2, 0.5]), step, 1e-7, 50) axs_gaps[1].semilogy(gaps_mirror, "-", marker=markers[i], color="purple", label=rf"$\gamma=${step}", ms=4, alpha=1-steps_norm[i]) axs_gaps[0].semilogy(gaps_pgd, "-", marker=markers[i], color="red", label=rf"$\gamma=${step}", ms=4, alpha=1-steps_norm[i]) axs_gaps[1].set_title("MD") axs_gaps[0].set_title("PGD") axs_gaps[0].legend(loc=1) axs_gaps[1].legend(loc=1) axs_gaps[0].set_xlabel("Itérations") axs_gaps[1].set_xlabel("Itérations") axs_gaps[0].set_ylabel(r"$Gap(x_n)$")fig_gaps.subplots_adjust(top=0.8)plt.suptitle(r"Evolution du saut de dualité ")plt.show()```#### b) Stabilité des itérésA partir de la liste des itérés, on peut monitorer la stabilité de la trajectoire des itérés.Un critère d'arrêt souvent utilisé en pratique est de s'arrêter à l'étape $k>1$ si:$$\|x_k - x_{k-1}\|_2\leq \epsilon \enspace.$$Pour $\gamma\in\{0.1, 0.5, 1\}$, visualiser l'évolution de ce critère pour PGD et MD. Commenter le résulat obtenu (on pourra s'aider en regardant le contenu de la liste itérés obtenus).```{python}#| code-fold: true#| fig-align: center#| label: fig-stab#| fig-cap: Évolution de la stabilités des itérés.steps = np.array([0.1, 0.5, 1])markers = ["o", "v", "s", "D", "X", "P"]steps_norm = np.linspace(0.1,0.95, len(steps))def stab_iteres(ll): norme = []for i, _ inenumerate(ll):if i ==0:continueelse: norme.append(np.linalg.norm(ll[i] - ll[i-1]))return normefig_gaps, axs_gaps = plt.subplots(1, 2, figsize=(8, 3), sharey=True)for i, step inenumerate(steps): ll_mirror, x_mirror, gaps_mirror = mirror_descent(grad_quad, np.array([0.2, 0.5]), step, 1e-7, 50) ll_pgd, x_pgd, gaps_pgd = projected_gd(grad_quad, proj_cons, np.array([0.2, 0.5]), step, 1e-7, 50) normes_md = stab_iteres(ll_mirror) normes_pgd = stab_iteres(ll_pgd) axs_gaps[1].semilogy(range(1, len(normes_md) +1), normes_md, "-", marker=markers[i], color="purple", label=rf"$\gamma=${step}", ms=4, alpha=1-steps_norm[i]) axs_gaps[0].semilogy(range(1, len(normes_pgd) +1), normes_pgd, "-", marker=markers[i], color="red", label=rf"$\gamma=${step}", ms=4, alpha=1-steps_norm[i]) axs_gaps[1].set_title("MD") axs_gaps[0].set_title("PGD") axs_gaps[0].legend(loc=1) axs_gaps[1].legend(loc=1) axs_gaps[0].set_xlabel("Itérations") axs_gaps[1].set_xlabel("Itérations") axs_gaps[0].set_ylabel(r"$||x_k - x_{k-1}||_2$")fig_gaps.subplots_adjust(top=0.8)plt.suptitle(r"Evolution de la stabilité des itérés")plt.show()```### Question 6 : Solution avec `scipy`?Utiliser la fonction `minimize` de `scipy` pour minimiser $f_{test}$ sur $[0,1]^2$.Comparer et commenter avec les résultats obtenus précédement.```{python}#| output: falsefrom scipy.optimize import minimizeminimize(f, np.array([0.2, 0.5]), method="TNC", bounds=((0,1),(0,1)))```### Question 7 : Résolution théorique du problèmeProuvez mathétatiquement quel est le minimum global de $f_{test}$ sur l'ensemble de contraintes $\mathcal{C}=[0,1]^2$.### Question 8 : Changeons l'espace de contraintesTrouver avec `scipy` la solution obtenue pour le problème suivant$$\argmin_{\mathcal{X}=(x_1,x_2)\in\mathbb{R}^2} f_{test}(x) \text{ sc } \begin{cases} y-x &< 0 \\ x^2 + 3 y^2 &\geq 0 \end{cases} \enspace.$$```{python}#| output: falseminimize(f, np.array([0.2, 0.5]), method="SLSQP", constraints = ( {'type': 'ineq', 'fun': lambda x: x[0] - x[1]}, {'type': 'ineq', 'fun': lambda x: x[0] **2+3* x[1] **2}))```# Optimisation non-convexeDans cette partie, nous utiliserons une variation de la fonction $f$ définie par$$f(x, y) = 0.9\sin(x) \times 0.9\cos(y) + 0.1 (x^2 + y^2) - 0.5$$Cette fonction n'est pas convexe, mais admet un minimum global !```{python}def tricky_func(x):"""Function with 2 local minima and one global minimum.""" x1, x2 = x[0], x[1]return0.9* np.sin(x1) *0.9* np.cos(x2) +0.1* (x1**2+ x2**2) -0.5def tricky_grad(x):"""Gradient of the function.""" x1, x2 = x[0], x[1] dfdx1 =0.9*0.9* np.cos(x1) *0.9* np.cos(x2) +0.2* x1 dfdx2 =-0.9* np.sin(x1) *0.9* np.sin(x2) +0.2* x2return np.array([dfdx1, dfdx2])def tricky_hess(x): x1, x2 = x[0], x[1] d2fdx1 =-0.9*0.9* np.sin(x1) *0.9* np.cos(x2) +0.2 d2fdx2 =-0.9* np.sin(x1) *0.9* np.cos(x2) +0.2 d2fdx1dx2 =-0.9*0.9* np.cos(x1) *0.9* np.sin(x2) d2fdx2dx1 = d2fdx1dx2return np.array([[d2fdx1, d2fdx1dx2], [d2fdx2dx1, d2fdx2]])```### Question 9 : Visualisation de la fonctionVisualiser la fonction en 3D ainsi que ses lignes de niveau.```{python}#| code-fold: true#| fig-align: center#| label: fig-Ackley#| fig-cap: Fonction et lignes de niveau.from mpl_toolkits.mplot3d import Axes3Dx = np.linspace(-5, 5, 100)X, Y = np.meshgrid(x, x)Z = np.zeros_like(X)for i inrange(len(x)):for j inrange(len(x)): Z[i,j] = tricky_func(np.array([X[i,j], Y[i,j]]))fig = plt.figure(figsize=(10, 3.6))ax1 = fig.add_subplot(121, projection='3d')ax1.plot_surface(X, Y, Z, cmap='RdBu_r')ax1.set_zlabel(r'$f(x)$')ax1.set_title('Visualisation de la fonction')# Plot the Ackley function with contourfax2 = fig.add_subplot(122)ax2.axes.set_aspect('equal')im = ax2.contourf(X, Y, Z, levels=30, cmap='RdBu_r')ax2.set_title('Visualisation des lignes de niveau')cbar = fig.colorbar(im, ax=ax2, orientation='vertical')plt.show()```### Question 10 : Minima et points d'initialisationPrendre $5$ points uniformément au hasard avec un générateur du module `np.random` fixé avec une seed fixée à $1$ dans le carré $[-3, 3]^2$ et éxécuter une descente de gradient sans contrainte avec un pas de $\gamma=0.1$, comparer les résultats obtenus. Commenter les minima obtenus et les problèmes rencontrés. Que faire en pratique quand on rencontre une fonction non convexe à minimiser (proposer deux solutions possibles en les commentant brièvement).```{python}#| code-fold: true#| fig-align: center#| label: fig-non-convex#| fig-cap: Initialization et convergence.rng = np.random.default_rng(1)points = rng.uniform(low=-3, high=3, size=(6, 2))ll_fin = []for i inrange(len(points)): ll_id, x_id, _ = projected_gd(tricky_grad, proj_id, points[i, :], 0.1, 1e-10, 500) ll_fin.append(ll_id)colors = ['b', 'g', 'r', 'c', 'm', 'y']x = np.linspace(-5, 5, 100)X, Y = np.meshgrid(x, x)fig, ax = plt.subplots(1,1)ax.contourf(X, Y, Z, levels=30, cmap='RdBu_r')ax.set_title('Lignes de niveau et itérations')for i inrange(len(points)): coords = np.array(ll_fin[i]) plt.plot(coords[:, 0], coords[:, 1], "o-", color=colors[i], label=f"Trajectoire #{i}")plt.legend()plt.show()```
