Introduction à la Recherche Opérationelle


M1 Statistique Science des Données (SSD)
Anatoli Juditsky
Franck Iutzeler
2017/2018

1. Problème de régression parcimonieuse et LASSO

Comme dans le TP sur le Dantzig Selector, on considère le modèle de régression normale $$ y=X\theta+\sigma\xi,\;\;\xi\sim \mathcal{N}(0, I_m),$$ où $X\in \mathbb{R}^{m\times n}$ et $y\in \mathbb{R}^m$ sont les observables, et $\theta\in \mathbb{R}^n$ est le paramètre inconnu. L'estimateur de $\theta$ par le LASSO ( least absolute shrinkage and selection operator cf. Tibshirani, Robert (1996). "Regression Shrinkage and Selection via the lasso". Journal of the Royal Statistical Society. Series B (methodological). Wiley. 58 (1): 267–88) peut être utilisé pour estimer $\theta$ dans le cas d'un modèle surparamétré, quand la dimension $n$ de $\theta$ est supérieure a la dimension de l'observation $y$. Dans ce cas l'estimateur ${\theta}_{L}$ s'écrit comme une solution du probleme d'optimisation

$$ \frac{1}{2} \| X\theta - y \|_2^2 + \lambda \|\theta\|_1 ,$$

où $\lambda>0$ est le paramètre de régularisation qui controle le degré de parcimonie de la solution.

Votre objectif dans cet exercice sera d'implementer l'estimateur ${\theta}_{L}$ en utilisant CVXPY et d'examiner comment la parcimonie des solutions évolue avec $\lambda$.

Question: Vérifier que le problème et les contraintes peuvent se formuler via des fonctions disponibles pour CVXPY.

2. Un exemple

In [ ]:
import cvxpy as cp
import numpy as np
In [ ]:
# Exemple jouet
n = 5
m = 3
sigma = 0.5

X = np.random.randn(m,n)
theta_true = np.array([1,0,0,0,5])
xi = np.random.randn(m)
y = np.dot(X,theta_true) + sigma*xi
In [ ]:
X
In [ ]:
theta_true
In [ ]:
y

Question: Trouver l'estimateur $\theta_{L}$ à partir de $X$ et $y$ par résolution du problème d'optimisation via CVXPY avec $\lambda$ fixé à 0.1.

In [ ]:
 

Question: Examiner la valeur de la solution pour $\lambda = 0$ puis $1$ puis $10$. Que remarquez vous sur la solution (coefficients proches de zéro, distance à la vraie solution)?

3. Fonction "Lasso"

Question: Écrivez une fonction MyLasso qui fait appel a CVXPy pour calculer l'estimation ${\theta}_{L}$.

Cette fonction doit sortir la liste avec les elements

  • coef, vecteur des coefficients de regression
  • resid, vecteur $y-X{\theta}_{L}$ de résidus
  • status, le champ status de la solution de CVXR.

L'appel à cette fonction devra être:

def MyLasso(X, y, lambda = 0.1)

  • X et y sont les observables
  • lambda est la valeur du paramètre de régularisation.
In [ ]:
 

a. Test quand on connait le vrai theta

Question: Testez votre fonction sur l' exemple ci-dessous

In [ ]:
m = 30
n = 20


X = np.random.randn(m,n)

theta_true = np.kron(np.array([0,0.5,1.0,4,20]),np.ones(int(n/5)))


xi = np.random.randn(m)
y = np.dot(X,theta_true) + sigma*xi

b. Variation du paramètre $\lambda$

Quand $\lambda=0$, peu de coefficients sont nul (ou presque, disons $<10^{-9}$); quand $\lambda \to \infty$, tous les coefficients deviennent nuls. Les solutions des valeurs intermédiaires forment le chemin de régularisation.

Question: Examinez les solutions du problème précédent pour $\lambda = 10^{-5}, 10^{-4}, ..,10^2$. Regardez notamment quelles valeurs sont nulles et comparer avec la vraie solution.

In [ ]:
 

c. Choix par validation croisée

Question:

  • Découpez $X$ et $y$ en une partie apprentissage ($X_{app}$ et $y_{app}$) et une partie test ($X_{test}$ et $y_{test}$).
  • Faites varier le paramètre de régularisation et pour chaque valeur:
    1. Calculez la solution $\theta_{L}$ du problème de lasso sur le jeu d'apprentissage ($\theta_L \in \arg\min_\theta \| X_{app}\theta - y_{app} \|_2^2 + \lambda_1 \|\theta \|_1 $ )
    1. Calculez l'erreur de cette solution $\| X_{test}\theta_L - y_{test} \|_2^2 $ sur le jeu de test
  • Déterminez une valeur pertinente du paramètre de régularisation par validation de l'error moyenne absolue sur la base de test
In [ ]: