Introduction à la Recherche Opérationelle


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



TP 4 - Classification


Problème de classification et sélection de variables

Nous considérons le modèle suivant de classification (prédiction de 0 et 1) $$ y=\mathrm{sign}(X\theta+\sigma\xi),\;\;\xi\sim \mathcal{N}(0, I_m),$$ où $X\in \mathbb{R}^{m\times n}$ et $y\in \{-1,1\}^m$ sont les observables, et $\theta\in \mathbb{R}^n$ est le paramètre inconnu. $\mathrm{sign}(x) = -1 $ si $x<0$ et $+1$ si $x\geq 0$.

Nous allons voir plusieurs méthodes pour prédire $y$ et les comparer. Les fonctions utilisées sont toutes des combinaisons des fonctions CVXPY disponibles ici: https://www.cvxpy.org/tutorial/functions/index.html#scalar-functions

Un exemple

In [1]:
import cvxpy as cp
import numpy as np
In [2]:
# Exemple jouet
n = 4
m = 6
sigma = 0.1

X = np.random.randn(m,n) # matrice de mélange -- connue

theta0 = np.random.randn(n)
non_null_coordinates = np.random.choice(n, int(0.8*n), replace=False)
theta_true = np.zeros(n)
theta_true[non_null_coordinates] = theta0[non_null_coordinates] # vrai theta -- inconnu: à estimer

xi = np.random.randn(m) # bruit 

y = np.sign(np.dot(X,theta_true) + sigma*xi) # observations -- connue
In [3]:
X
Out[3]:
array([[-0.07962712,  0.52175377, -0.34780377,  0.16200015],
       [-0.54117825,  0.99208539, -1.57374623,  0.70350946],
       [ 0.34264479, -0.52427128, -1.4801209 , -0.09511679],
       [ 0.55082506,  1.77925949, -0.72487898, -1.12467504],
       [ 0.36551622,  1.57705249,  1.29123143, -1.5245815 ],
       [-0.25885727, -0.12820224,  0.52597023,  0.07359679]])
In [4]:
theta_true
Out[4]:
array([-0.02985106,  0.        , -0.19189802,  0.79324613])
In [5]:
y
Out[5]:
array([ 1.,  1.,  1., -1., -1.,  1.])

Estimateur A

Le premier estimateur que nous allons implémenter est le lasso:

\begin{equation} \tag{A} \theta_A \in \arg\min_\theta \left\{ \frac{1}{2} \| X\theta - y \|_2^2 + \lambda \|\theta\|_1 \right\}, \end{equation}

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

Question 1: Écrivez une fonction EstimateurA qui retourne ${\theta}_{A}$.

  • Cette fonction qui prend en entrée:
    • X
    • y
    • lambda
  • et sort les elements
    • coef, vecteur des coefficients de regression
    • status, le champ status de la solution du problème d'optimisation.

Testez cette fonction sur l'exemple du début.

In [ ]:

Estimateur B

Prenons un estimateur proche du sélecteur de Dantzig (TP2)

\begin{equation} \tag{B} \theta_B \in \arg\min_\theta \left\{ \|\theta\|_1 \text{ sous contrainte } \| X\theta - y \|_2^2 \leq \lambda \right\}, \end{equation}

où $\lambda>0$ est l'erreur autorisée qui restreint les solutions à être optimale à $\lambda$ près.

Question 2: Écrivez une fonction EstimateurB qui retourne ${\theta}_{B}$.

  • Cette fonction qui prend en entrée:
    • X
    • y
    • lambda
  • et sort les elements
    • coef, vecteur des coefficients de regression
    • status, le champ status de la solution du problème d'optimisation.

Testez cette fonction sur l'exemple du début.

In [ ]:

Question 3: Testez l'estimateur B pour différentes valeurs de lambda. Le problème d'optimisation a-t-il toujours une solution? Expliquez.

In [ ]:

Estimateur C

Au lieu de contraindre l'erreur, contraignons les valeurs du paramètre.

\begin{equation} \tag{C} \theta_C \in \arg\min_\theta \left\{ \| X\theta - y \|_2^2 \text{ sous contrainte } |\theta_i| \leq \lambda \text{ pour tout } i \right\}, \end{equation}

où $\lambda>0$ est la norme maximale.

Question 4: Écrivez une fonction EstimateurC qui retourne ${\theta}_{C}$.

  • Cette fonction qui prend en entrée:
    • X
    • y
    • lambda
  • et sort les elements
    • coef, vecteur des coefficients de regression
    • status, le champ status de la solution du problème d'optimisation.
    • dual, la variable duale optimale

Testez cette fonction sur l'exemple du début.

In [ ]:

Question 5: Faites varier lambda et utilisez la valeur de la variable duale pour déterminer si la contrainte est active ou inactive.

In [ ]:

Question 6: La solution $\theta_C$ est-elle toujours différente de la solution du problème des moindres carrés $\theta_{MoindresCarrés} \in \arg\min_\theta \left\{ \| X\theta - y \|_2^2 \right\}$ ?

In [ ]:

Comparison des estimateurs

Comparons maintenant les estimateurs ci-dessus à partir de l'exemple suivant:

In [6]:
# Exemple plus grand
n = 60
m = 200
sigma = 0.1

X = np.random.randn(m,n) # matrice de mélange -- connue

theta0 = np.random.randn(n)
non_null_coordinates = np.random.choice(n, int(0.8*n), replace=False)
theta_true = np.zeros(n)
theta_true[non_null_coordinates] = theta0[non_null_coordinates] # vrai theta -- inconnu: à estimer

xi = np.random.randn(m) # bruit 

y = np.sign(np.dot(X,theta_true) + sigma*xi) # observations -- connue

Divisons les données à moitié en jeu d'entrainement et moitié en jeu de test.

In [7]:
X_train = X[:int(m/2),:]
y_train = y[:int(m/2)]

X_test = X[int(m/2):,:]
y_test = y[int(m/2):]

Question 7: Pour chacun des estimateurs, trouvez la valeur de $\lambda$ pour laquelle l'estimateur (calculé sur les données d'entrainement) minimise l'erreur de classification (sur les données de test).

Ex: Pour l'estimateur A, on calcule theta_A, = EstimateurA(X_train, y_train, lam = 1e-5 ) puis la prédiction $$ y_{A,pred} = \text{sign}(X_{test}\theta_A). $$ Puis l'erreur de classification s'écrit comme le pourcentage de coordonnées $i$ pour lesquelles $ y_{A,pred}(i) \neq y_{test}(i)$.

In [ ]:

Question 8: Comparez les quatre estimateurs. (Cette question est volontairement ouverte, plusieurs critères peuvent être utilisés, des propositions d'améliorations peuvent etre formulées.)

In [ ]: