Introduction à la Recherche Opérationelle


M1 Statistique Science des Données (SSD)
Anatoli Juditsky
Franck Iutzeler
2017/2018

CVXR est un package R permettant la résolution de problèmes d'optimisation une fois ceux-ci formulés correctement. Il est disponible sur le serveur CRAN et s'installe donc simplement par

    install.packages("CVXR")

Le site de CVXR https://cvxr.rbind.io/ fournit de nombreuses ressources, notamment

mais plus important la liste des fonctions disponibles.

In [ ]:
library(CVXR)

1. Un premier exemple

Considérons le problème de régression linéaire $$ Y=X\beta^\star+\epsilon$$ où $Y$ est un vecteur de taille 20, $X$ est une matrice de taille 20x10 et $\beta^\star =[-4,..,-1,0,1,..,5]$ est un vecteur de taille 10; finalement, epsilon est un vecteur aléatoire Gaussien de taille 20 et de variance 0.01.

In [2]:
n <- 20
p <- 10
beta_star <- -4:5   # beta is just -4 through 5.

X <- matrix(rnorm(n * p), nrow=n)
Y <- X %*% beta_star + 0.1*rnorm(n)

L'estimateur du maximum de vraisemblance de $\hat{\beta}$ sachant $Y$ et $X$ est le vecteur minimisant en $\beta$ la fonction objective $$ L(\beta) = \| X \beta - Y \|_2^2 $$

 Resolution avec CVXR

  • Étape 1. Définir la variable selon laquelle minimiser
In [3]:
beta <- Variable(p) # p est la taille de la variable, définie ci-dessus
  • Étape 2. Définir la fonction objectif
In [4]:
L <- sum(( X %*% beta  -  Y)^2) # L est la fonction ci-dessus 

objective <- Minimize(L)

Ici, la variable L est bien une fonction pour CVXR comme son expression fait intervenir un objet beta qui est une variable CVXR.

L'objectif est ainsi bien de minimiser L (sous entendu par rapport à la variable beta ).

  • Étape 3. Créer le problème à résoudre
In [5]:
problem <- Problem(objective)
  • Étape 4. Le résoudre!
In [6]:
result <- solve(problem)
  • Étape 5. Analyser le résultat

La variable de sortie du solver possède (entre autres) les attributs suivants:

$status : est-ce que la solution trouvée est optimale

$value : valeur optimale de L

$solve_time : le temps mis à trouver une solution

$getValue : la fonction permettant d'avoir la valeur optimale (de la variable précisée)

In [7]:
result$status
'optimal'
In [8]:
result$solve_time
0.000335341
In [9]:
betaHat <- result$getValue(beta)
print(betaHat)
print(beta_star)
             [,1]
 [1,] -3.96121648
 [2,] -2.98943014
 [3,] -2.00700971
 [4,] -1.04417007
 [5,]  0.01668267
 [6,]  1.00525186
 [7,]  2.02863310
 [8,]  2.97801077
 [9,]  3.96458272
[10,]  5.03132346
 [1] -4 -3 -2 -1  0  1  2  3  4  5

Conclusion

On a donc résolu le problème de régression linéaire assez facilement, comme on aurait pu le faire avec le package lm. Cependant, comme nous alons le voir maintenant, CVXR permet aisément de modifier le problème, en ajoutant par exemple des contraintes.

2. Un problème avec contraintes

Supposons que nous voulons maintenant résoudre le même problème de minimiser $$ L(\beta) = \| X \beta - Y \|_2^2 $$ mais avec la contrainte supplémentaire que les composantes de $\beta$ doivent être positives et que leur somme doit être égale à 10.

Suivons les même étapes que précedemment.

  • Étape 1. Définir la variable selon laquelle minimiser
In [10]:
beta <- Variable(p) # p est la taille de la variable, définie ci-dessus
  • Étape 2. Définir la fonction objectif et les contraintes
In [11]:
L <- sum(( X %*% beta  -  Y)^2) # L est la fonction ci-dessus 

objective <- Minimize(L)
In [12]:
contrainte1 <- beta>=0

un <- rep(1,p)
contrainte2 <- un%*%beta == 10
  • Étape 3. Créer le problème à résoudre avec les contraintes
In [13]:
problem <- Problem(objective , constraints = list(contrainte1, contrainte2))
  • Étape 4. Le résoudre!
In [14]:
result <- solve(problem)
  • Étape 5. Analyser le résultat
In [15]:
result$status
'optimal'
In [16]:
betaHat2 <- result$getValue(beta)
betaHat2
sum(betaHat2)
-5.930560e-09
-8.454154e-09
1.722454e-08
-5.130932e-09
3.906343e-08
2.408008e-08
2.581076e+00
3.752495e-07
3.081011e+00
4.337912e+00
9.99999999985069

Le vecteur obtenu vérifie bien les contraintes ajoutées.

(il est également possible de récupérer les variables duales associées aux contraintes)

In [17]:
result$getDualValue(contrainte1)
1.564415e+02
1.537000e+02
3.431968e+01
1.164025e+02
2.003450e+01
3.239110e+01
4.186425e-07
5.904233e+00
3.525306e-07
2.530621e-07
In [18]:
result$getDualValue(contrainte2)
17.486351481256