Descente de gradient projeté

Auteur·rice

Joseph Salmon, Tanguy Lefort

Objectifs de ce TP
  • Comment optimiser une fonction sous contraintes? Application au problème des moindres carrés linéaires avec données réelles. Application du théorème de projection sur un ensemble convexe fermé. Proposer un modèle pour des données réelles.
  • Implémenter une projection. Combiner deux algorithmes (projection et descente de gradient). Charger des données réelles dans Python. Visualiser des données réelles. Utiliser la librairie scikit-learn.

Dans ce TP, nous allons utiliser la méthode de descente de gradient pour trouver un minimiseur d’une fonction f:\mathbb{R}^d\rightarrow\mathbb{R} pour d\geq 2 en posant des contraintes sur f ou sur l’ensemble \mathcal{C} des solutions du problème.

Algorithme de la descente de gradient projeté

Dans cette partie, on suppose que l’on a une fonction f, et on va introduire la descente de gradient projeté (🇬🇧 projected gradient descent) avec l’ensemble de contraintes convexe fermé \mathcal{C}. On note \Pi_{\mathcal{C}} la projection convexe, i.e. LA1 fonction convexe qui envoie un élément dans l’ensemble convexe fermé \mathcal{C}.

    \begin{algorithm}
    \caption{Descente de gradient projeté (à pas fixe)}
    \begin{algorithmic}
    \INPUT $\nabla f, \Pi_{\mathcal{C}}, x_{init}, \gamma,$ maxiter, $\epsilon$
    \STATE $x = x_{init}$
    \COMMENT{Initialization}
    \FOR{$i=1,\dots,$maxiter}
        \STATE $g\leftarrow \nabla f(x)$
        \COMMENT{Evaluation du gradient en l'itéré courant}
    \IF{$\|g\|^2 \leq \epsilon^2$ }
        \COMMENT{Caractérisation d'un extremum}
        \STATE Break
    \ELSE
    \STATE $x \leftarrow \Pi_{\mathcal{C}}\big(x - \gamma g\big)$
    \COMMENT{Pas dans la direction opposée au gradient}
    \ENDIF
    \ENDFOR
    \OUTPUT $x$
    \end{algorithmic}
    \end{algorithm}

Dans la suite de cette partie, on travaillera avec la fonction f_{test} définie ainsi:

\begin{align*} f_{test}:\quad \mathbb{R}^2 &\to \mathbb{R}^2\\ \begin{pmatrix} x_1 \\ x_2 \\\end{pmatrix} & \mapsto (x_1 - 1)^2 + 3(x_2 + 1)^2 \end{align*}

Question 1: Minimiser avec ou sans contraintes

Quel est le minimum global de f_{test} ? (on pourra utiliser la descente de gradient du TP précédent). Utiliser la méthode minimize de la librairie scipy pour trouver le minimum de f_{test} sur \mathbb{R}^2_+. Pour cela, utiliser l’argument bounds et une méthode d’optimisation qui supporte cette contrainte, par exemple method="TNC". Voir la documentation si besoin.

Question 2: Descente de gradient projeté

Implémenter la méthode de descente du gradient projeté à pas fixe en prenant comme argument une fonction de projection \Pi_{\mathcal{C}}:\mathbb{R}^d\to \mathbb{R}^d quelconque.

Question 3: Application à la contrainte de positivité

Avec l’algorithme précédent, résoudre numériquement : \begin{aligned} \min_{x\in\mathbb{R}^2}\ & f_{\text{test}}(x) \\ & \text{ s.c. } \begin{cases} x_1 \geq 0 \\ x_2 \geq 0 \end{cases} \end{aligned}

Note

Cette formulation induit la positivité des paramètres. Elle est notamment utilisée en physique quand il y a une interprétation directe comme la vitesse ou la masse d’un objet.

Question 4: Application aux projections sur les boules

Avec l’algorithme précédent, résoudre numériquement : \begin{aligned} \min_{x\in\mathbb{R}^2}\ & f_{\text{test}}(x) \\ & \text{ s.c. } \|x\|_2^2\leq 1 \end{aligned} Adapter l’exemple ci-dessus pour résoudre avec la contrainte \|x\|^2\leq 1/2. Vous pouvez dans un premier temps tester votre fonction de projection sur quelques points dont les projetés sont faciles à calculer à la main.

Pour aller plus loin

En s’aidant de l’exemple dans la documentation de la fonction minimize de scipy, résoudre le problème suivant : \begin{aligned} \min_{x\in\mathbb{R}^2}\ f_{\text{test}}(x) \text{ s.c. } \begin{cases} x_1 - 3 x_2 \geq -2 \\ x_1 + x_2 \geq 0 \end{cases} \enspace. \end{aligned} Ce problème est avec des contraintes d’inégalités linéaires. Pour le résoudre à la main, il serait possible d’utiliser la méthode des multiplicateurs de Lagrange.

Problème des moindres carrés linéaires

Dans cette partie, le but est de trouver la fonction linéaire qui approche le mieux un nuage de points (au sens de la norme euclidienne). C’est-à-dire, pour un ensemble de n points de \mathbb{R^2} (x_i, y_i)_{i=1}^n nous cherchons à résoudre le problème convexe suivant:

\hat w = \begin{pmatrix} \hat w_0 \\ \hat w_1\end{pmatrix}\in \argmin_{w\in\mathbb{R}^2} \underbrace{\frac{1}{2n} \sum_{i=1}^n (y_i - (w_0 + w_1 x_i))^2}_{f(w)} \enspace.

Une interpretation est donnée ci-dessous: Visualisation de la méthode des moindres carrés

Question 1: Moindres carrés

Pour résoudre le problème des moindres carrés linéaires, utilisons la descente de gradient pour obtenir \hat w. Pour cela, on notera que l’on peut écrire f(w)=\frac{1}{2n}\|y-Xw\|_2^2 avec:

y=\begin{pmatrix} y_1 \\ \vdots \\ y_n\end{pmatrix}\in\mathbb{R}^n,\quad X=\begin{pmatrix} 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n\end{pmatrix}\in\mathbb{R}^{n\times 2} \enspace.

Il faut donc calculer (mathématiquement) le gradient \nabla f(w) puis appliquer la descente de gradient pour résoudre le problème.

Implémenter une fonction gradient qui prend en entrée y,X,w et renvoie \nabla f(w) (ne pas utiliser explicitement le fait que X ait seulement 2 colonnes pour pouvoir généraliser par la suite à plus de paramètres).

Question 2: Visualiser les données

En utilisant le jeu de données Iowa_Liquor_tp.csv et le code ci-dessous pour charger un fichier .csv, tracer le nuage de points (x_i, y_i)_{i=1}^{n}.

import pandas as pd

# Whole dataset: https://data.iowa.gov/api/views/m3tr-qhgy/rows.csv?accessType=DOWNLOAD
# dataset info: https://data.iowa.gov/Sales-Distribution/Iowa-Liquor-Sales/m3tr-qhgy

# Warning for MAC OS users: https://stackoverflow.com/questions/40684543/how-to-make-python-use-ca-certificates-from-mac-os-truststore

url = "https://tanglef.github.io/_teachings/data/HAX606X_optim/datasets/Iowa_Liquor_tp.csv"
data = pd.read_csv(url)
data.rename(columns={"Unnamed: 0": "Date"}, inplace=True)
print(data)
           Date  Bottles Sold
0    2012-01-01       1302231
1    2012-02-01       1457824
2    2012-03-01       1458665
3    2012-04-01       1588914
4    2012-05-01       1782175
..          ...           ...
114  2021-07-01       2558550
115  2021-08-01       2646637
116  2021-09-01       2589685
117  2021-10-01       2709511
118  2021-11-01       2619605

[119 rows x 2 columns]

Pour accéder aux valeurs d’une colonne, on fera comme dans l’exemple ci-dessous:

data["Bottles Sold"]
0      1302231
1      1457824
2      1458665
3      1588914
4      1782175
        ...   
114    2558550
115    2646637
116    2589685
117    2709511
118    2619605
Name: Bottles Sold, Length: 119, dtype: int64

Données brutes et régressions

from sklearn import linear_model
import numpy as np


y = data["Bottles Sold"].values.reshape(-1, 1)
x = np.arange(len(y)).reshape(-1, 1)

lm = linear_model.LinearRegression()
lm.fit(x, y)
predictions = lm.predict(x)


x2 = np.sin(2 * np.pi / 6 * x)
x3 = np.cos(2 * np.pi / 6 * x)
X = np.hstack((np.ones_like(x), x, x2, x3))

lm_period = linear_model.LinearRegression(fit_intercept=False)
lm_period.fit(X, y)
predictions_period = lm_period.predict(X)

On observe qu’il y a une périodicité annuelle dans les données. Implémenter alors la même méthode pour en tenir compte. Pour cela, résoudre le problème suivant:

\begin{pmatrix} \hat{w}_0 \\ \hat{w}_1 \\ \hat{w}_2 \\ \hat{w}_3 \end{pmatrix} \in \argmin_{w\in\mathbb{R}^4} \underbrace{ \sum_{i=1}^n \Big[y_i - \left(w_0 + w_1 \cdot x_i + w_2 \cdot \cos\big(\tfrac{2 \pi x_i}{6}\big) + w_3 \cdot \sin\big(\tfrac{2 \pi x_i}{6}\big) \right) \Big]^2}_{f(w_0, w_1, w_2, w_3)}

Vérifiez que votre solution coïncide avec la solution proposée par la librairie scikit-learn avec le code suivant:

what = lm.coef_
w0hat = lm.intercept_  # ordonnée à l'origine
print(what, w0hat)
[[9226.57063096]] [1464946.88739496]
Pour aller plus loin

La méthode pour résoudre le problème des moindres carrés proposée est celle de scikit-learn. Cette librairie est extrêmement utilisée en apprentissage statistique. Explorez la galerie de la documentation et découvrez plein de nouvelles méthodes et d’applications!