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 [1]:
require(CVXR)
Loading required package: CVXR

Attaching package: ‘CVXR’

The following object is masked from ‘package:stats’:

    power

In [2]:
# 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 [3]:
X
theta_true
y
0.1793983 -1.0645670 -0.46074551-1.1719733 -0.3459315
0.4851292 0.5802169 -0.04891822 0.0672531 0.8001704
-1.4484111 0.6131263 -0.17330022 0.4770147 0.4234673
  1. 1
  2. 0
  3. 0
  4. 0
  5. 5
-1.4490105
4.5868486
0.6208041

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.

In [4]:
theta <- Variable(n)
In [5]:
F <- cvxr_norm(theta, 1)
objective <- Minimize(F)
In [6]:
kappa = 0.2
contrainte <- cvxr_norm( t(X)%*%(X%*%theta - y) , "inf") <= kappa*sigma
In [7]:
problem <- Problem(objective , constraints = list(contrainte))
In [8]:
result <- solve(problem)
In [9]:
result$status
theta_ds=result$getValue(theta)
'optimal'
In [10]:
cbind( theta_ds , theta_true)
theta_true
1.026016e+001
7.677025e-120
-1.218286e-110
-8.879506e-020
5.076118e+005

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).$
In [11]:
DSelect <- function(X, y, sigma = 1, c = 1, rtol = 1e-6, verb = 5){

    d = dim(X)
    m = d[1]
    n = d[2]
     
    alpha = 0.05
    kappa = c*qnorm(1-alpha/m, 0, 1)

    theta <- Variable(n) 
    # Function
    F <- cvxr_norm(theta, 1)
    objective <- Minimize(F)

    # Constraints
    contrainte <- cvxr_norm( t(X)%*%(X%*%theta - y) , "inf") <= kappa*sigma
    
    problem <- Problem(objective , constraints = list(contrainte))
    
    result <- solve(problem)

    theta_ds=result$getValue(theta)
    
    DS = list()
    DS$coef = theta_ds
    DS$resid =  y - X%*% theta_ds
    DS$status = result$status
    return(DS)
}

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

a. Test quand on connait le vrai theta

In [12]:
n = 5
m = 5
X = replicate(n, rnorm(m))
theta_true = rnorm(n)
sigma = 0.005
y = X%*%theta_true + sigma*rnorm(m)
In [13]:
f = DSelect(X, y, sigma=sigma)
str(f)
List of 3
 $ coef  : num [1:5, 1] 1.42 3.09e-12 -1.96e-01 7.17e-01 -4.34e-01
 $ resid : num [1:5, 1] -0.013792 -0.006892 0.000163 0.007348 0.008817
 $ status: chr "optimal"
In [14]:
comp = cbind(as.matrix(f$coef), as.matrix(theta_true))
comp
1.422329e+00 1.3561656
3.085377e-12 0.1035487
-1.963468e-01-0.2585645
7.172774e-01 0.7746473
-4.338057e-01-0.3714834

b. Example du papier de Candes/Tao

In [15]:
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 [16]:
f <- DSelect(X, y, sigma= sigma)
In [17]:
S_set

plot(beta,col="red")
points(f$coef,col="blue")
  1. 67
  2. 55
  3. 51
  4. 173
  5. 219
  6. 61
  7. 224
  8. 14