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 [1]:
import cvxpy as cp
import numpy as np
import pandas as pd
In [2]:
# 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 [3]:
X
Out[3]:
array([[-1.16395435,  0.94891256,  0.14510104, -1.02031382,  1.13207932],
       [ 0.49051027,  0.20326417, -0.45235022,  0.16953611,  0.35772628],
       [ 0.12238295,  0.79245213, -1.26381771,  1.54739771,  1.82549954]])
In [4]:
y
Out[4]:
array([4.50486661, 2.15739625, 9.36706893])
In [5]:
theta_true
Out[5]:
array([1, 0, 0, 0, 5])

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 [6]:
theta = cp.Variable(n) 
In [7]:
F  = cp.norm1(theta)
objective = cp.Minimize(F)
In [8]:
kappa = 0.2
contrainte  = cp.norm_inf( np.dot(X.T,X)*theta - np.dot(X.T,y) ) <= kappa*sigma
In [9]:
problem = cp.Problem(objective , constraints = [contrainte])
In [10]:
problem.solve(verbose=True,solver="ECOS")
ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -2.000e-02  +1e+02  3e-01  9e-01  1e+00  4e+00    ---    ---    1  1  - |  -  - 
 1  +5.408e+00  +5.902e+00  +3e+01  7e-02  2e-01  7e-01  1e+00  0.8146  1e-01   0  0  0 |  0  0
 2  +5.045e+00  +5.135e+00  +6e+00  1e-02  3e-02  1e-01  3e-01  0.8381  6e-02   0  0  0 |  0  0
 3  +5.301e+00  +5.326e+00  +1e+00  2e-03  7e-03  3e-02  7e-02  0.7866  4e-02   0  0  0 |  0  0
 4  +5.475e+00  +5.502e+00  +8e-01  3e-03  5e-03  3e-02  4e-02  0.6882  4e-01   0  0  0 |  0  0
 5  +5.518e+00  +5.518e+00  +1e-01  3e-04  7e-04  2e-03  6e-03  0.9310  1e-01   0  0  0 |  0  0
 6  +5.531e+00  +5.531e+00  +2e-02  4e-05  1e-04  2e-04  1e-03  0.8357  4e-03   1  0  0 |  0  0
 7  +5.534e+00  +5.534e+00  +3e-03  5e-06  1e-05  4e-05  1e-04  0.9691  9e-02   1  0  0 |  0  0
 8  +5.534e+00  +5.534e+00  +3e-05  6e-08  1e-07  5e-07  1e-06  0.9890  1e-04   1  0  0 |  0  0
 9  +5.534e+00  +5.534e+00  +3e-07  6e-10  2e-09  5e-09  2e-08  0.9890  1e-04   1  0  0 |  0  0
10  +5.534e+00  +5.534e+00  +4e-09  7e-12  2e-11  6e-11  2e-10  0.9890  1e-04   1  0  0 |  0  0

OPTIMAL (within feastol=1.7e-11, reltol=6.4e-10, abstol=3.5e-09).
Runtime: 0.000680 seconds.

Out[10]:
5.534047876210565
In [11]:
problem.status
Out[11]:
'optimal'
In [12]:
theta_ds= theta.value
theta_ds
Out[12]:
array([ 1.11376725e-10,  5.88467640e-11, -1.24499619e+00,  8.52837887e-02,
        4.20376789e+00])
In [13]:
print("True \t Dantzig selector")
for i in range(len(theta_true)):
    print('{} \t {:6.3f}'.format(theta_true[i],theta_ds[i]))
True 	 Dantzig selector
1 	  0.000
0 	  0.000
0 	 -1.245
0 	  0.085
5 	  4.204

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 [14]:
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)

    theta = cp.Variable(n) 
    # Function
    F = cp.norm1(theta)
    objective = cp.Minimize(F)

    # Constraints
    contrainte  = cp.norm_inf( np.dot(X.T,X)*theta - np.dot(X.T,y) ) <= kappa*sigma
    
    problem = cp.Problem(objective , constraints = [contrainte])
    
    problem.solve(verbose=verb,solver="ECOS")

    theta_ds = theta.value
    
    residual =  y - np.dot(X,theta_ds)
    sol_status = problem.status
    
    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 [15]:
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 [16]:
f = DSelect(X, y, sigma=sigma)
f
Out[16]:
(array([-0.24718512, -1.22342349, -0.37820185,  0.26991672, -0.09681672]),
 array([ 0.00272682,  0.02139071, -0.00994895, -0.00415782, -0.00527264]),
 'optimal')
In [17]:
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]))
True 	 Dantzig selector
-0.231 	 -0.247
-1.254 	 -1.223
-0.384 	 -0.378
 0.265 	  0.270
-0.108 	 -0.097

b. Example du papier de Candes/Tao

In [18]:
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 [19]:
f = DSelect(X, y, sigma= sigma)
In [20]:
S_set
Out[20]:
[134, 133, 154, 170, 115, 211, 28, 120]
In [22]:
import matplotlib.pyplot as plt

plt.plot(beta,'ro')
plt.plot(f[0],'b*')
Out[22]:
[<matplotlib.lines.Line2D at 0x7fad2463e550>]