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

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.

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.

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]