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

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

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_{DS}$ à partir de $X$ et $y$ par résolution du problème d'optimisation via CVXR avec $\kappa$ fixé à 0.2.

3. Fonction "Dantzig Selector"

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

Cette fonction doit sortir la liste avec les elements

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

L'appel à cette fonction devra être:

DSelect <- function(X, y, sigma = 1, c = 1, rtol = 1e-6, verb = 5)

  • 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).$

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

a. Test quand on connait le vrai theta

In [ ]:
n = 5
m = 5
X = replicate(n, rnorm(m))
theta_true = rnorm(n)
sigma = 0.005
y = X%*%theta_true + sigma*rnorm(m)
In [ ]:
f = DSelect(X, y, sigma=sigma)
str(f)
In [ ]:
comp = cbind(as.matrix(f$coef), as.matrix(theta_true))
comp

b. Example du papier de Candes/Tao

In [ ]:
n = 256
m = 72
S = 8

S_set = sample.int(n, size = S, replace = FALSE)

X = replicate(n, rnorm(m))

beta = as.matrix(rep(0,n))
for (i in S_set){
 beta[i] = rnorm(1)
    }

sigma = 1/3.0 * sqrt(S/m)
y = X%*%beta + sigma*rnorm(m)
In [ ]:
f <- DSelect(X, y, sigma= sigma)
In [ ]:
S_set

plot(beta,col="red")
points(f$coef,col="blue")