Introduction à la Recherche Opérationelle


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

1. Problème de régression parcimonieuse et 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$ de ``Dantzig Selector'' (cf. Candes, E., Tao, T. (2007). The Dantzig selector: Statistical estimation when $p$ is much larger than $n$. The Annals of Statistics, 2313-2351) 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}_{DS}$ s'écrit comme une solution du probleme d'optimisation

$${\theta}_{DS}\in \arg\min_{\theta\in \mathbb{R}^n} \left\{\|\theta\|_1,\;\mbox{sous contrainte}\;\|X^T(X\theta-y)\|_\infty\leq \kappa\sigma\right\},$$

où $\kappa>0$ est un hyper-paramètre. La valeur de $\kappa$, préconisée dans la literature, est $c q_{\mathcal{N}}\left(1-{\alpha\over 2m}\right)$, où $\alpha\in (0,1)$ est le niveau de risque choisi (e.g., $\alpha=.05$, etc) et $q_\mathcal{N}(p)$ est la $p$-quantile de la normale standardisée, et $c=\max_j\|[X]_j\|_2=\max_j\sqrt{[X]_j^T[X]_j}$ est la norme maximale de colonne de la matrice $X$.

Votre objectif dans cet exercice sera d'implementer l'estimateur ${\theta}_{DS}$ en utilisant CVXPY.

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
import pandas as pd
In [ ]:
# Exemple jouet
n = 5
m = 3
sigma = 0.1

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 [ ]:
y
In [ ]:
theta_true

Question: Trouver l'estimateur $\theta_{DS}$ à partir de $X$ et $y$ par résolution du problème d'optimisation via CVXPY avec $\kappa$ fixé à 0.2 et utiliser le solver ECOS.

In [ ]:
 

3. Fonction "Dantzig Selector"

Question: Écrivez une fonction DSelect qui fait appel a CVXPY pour calculer l'estimation ${\theta}_{DS}$.

Cette fonction doit sortir un tuple avec les elements

  • coef, vecteur des coefficients de regression
  • resid, vecteur $y-X{\theta}_{DS}$ de résidus
  • status, le statut de sortie du solver

L'appel à cette fonction devra être:

DSelect(X, y, sigma = 1, c = 1, verb = False)

  • X et y sont les observables
  • sigma est une estimation de $\sigma$
  • c est le paramètre réel tel que la valeur de $\kappa$ dans {DS} est $ \kappa=c\,q_{\mathcal{N}}\left(1-{\alpha\over 2m}\right).$
In [ ]:
import scipy.stats

def DSelect(X, y, sigma = 1, c = 1, verb = False):
    

    m,n = X.shape
     
    alpha = 0.05
    kappa = c*scipy.stats.norm.ppf(1-alpha/m)

    #### TODO
    
    return theta_ds,residual,sol_status

Question: Testez votre fonction sur les deux exemples ci-dessous

a. Test quand on connait le vrai theta

In [ ]:
n = 5
m = 5
X = np.random.randn(m,n)
theta_true = np.random.randn(n)
sigma = 0.005

y = np.dot(X,theta_true) + sigma*np.random.randn(m)
In [ ]:
f = DSelect(X, y, sigma=sigma)
f
In [ ]:
print("True \t Dantzig selector")
for i in range(len(theta_true)):
    print('{:6.3f} \t {:6.3f}'.format(theta_true[i],f[0][i]))

b. Example du papier de Candes/Tao

In [7]:
import random

n = 256
m = 72
S = 8

S_set = random.sample(range(n),k=S)

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

beta = np.zeros(n)
beta[S_set] = np.random.randn(S)

sigma = 1/3.0*np.sqrt(S/m)
y = np.dot(X,beta) + sigma*np.random.randn(m)
In [ ]:
f = DSelect(X, y, sigma= sigma)
In [ ]:
S_set
In [ ]:
import matplotlib.pyplot as plt

plt.plot(beta,'ro')
plt.plot(f[0],'b*')