Introduction à la Recherche Opérationelle


M2 Statistique Science des Données (SSD)
Anatoli Juditsky
Franck Iutzeler
Version R2019

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 CVXR 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 CVXR.

2. Un exemple

In [ ]:
require(CVXR)
In [ ]:
# Exemple jouet
n = 5
m = 3
sigma = 0.5

X = replicate(n, rnorm(m))
theta_true = c(1,0,0,0,5)
xi = rnorm(m)
y = X%*%theta_true + sigma*xi
In [ ]:
X
theta_true
y

Question: Trouver l'estimateur $\theta_{L}$ à partir de $X$ et $y$ par résolution du problème d'optimisation via CVXR 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 CVXR 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:

MyLasso <- function(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 = replicate(n, rnorm(m))
theta_true = c( rep(0,floor(n/5)) , rep(0.5,floor(n/5)) , rep(1.0,floor(n/5)), rep(4,floor(n/5)), rep(20,floor(n/5)) )
sigma = 0.1
y = X%*%theta_true + sigma*rnorm(m)
In [ ]:
f = MyLasso(X, y)
str(f)
In [ ]:
comp = cbind(as.matrix(f$coef), as.matrix(theta_true))
comp

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 [ ]:

4. Le Lasso en apprentissage

Il s'agit dans cette partie de prédire les facteurs socio-économiques à partir d'une base de données issue des données du FBI de 1995 pour 1994 quartiers américains.

  • Les 5 premières colonnes ne sont pas prédictives:

      state: US state (by number) - not predictive
      county: numeric code for county - not predictive, and many missing values (numeric)
      community: numeric code for community - not predictive and many missing values (numeric)
      communityname: community name - not predictive - for information only (string)
      fold: fold number for non-random 10 fold 
  • Les 122 suivantes sont prédictives et numériques

      données socio-économiques (voir 'data/communities.name')
  • La dernière colonne contient la valeur à prédire

      ViolentCrimesPerPop: total number of violent crimes per 100K popuation (numeric - decimal) GOAL attribute (to be predicted)
  • Il y a des données manquantes

      Symbolisées par ? , vous pouvez les remplacer par 0
      vous devez faire attention que votre matrice finale est bien numérique (mode(X) = "numeric") pour y remédier) 

Question:

  • Lisez le dataset, découpez le en base d'apprentissage et de test.
  • Faites varier le paramètre de régularisation pour sélectionner les valeurs les plus importantes et les moins importantes selon le crière de mise à zéro par lasso par validation sur la base 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 [ ]:
crime_data = read.csv('./data/communities.data',sep=',', header=FALSE)
m = dim(crime_data)[1]
In [ ]: