Introduction à la Recherche Opérationelle


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

CVXPY est un package Python permettant la résolution de problèmes d'optimisation une fois ceux-ci formulés correctement. Il s'installe par pip

    pip install [--user] cvxpy

Le site de CVXPY https://www.cvxpy.org/ fournit de nombreuses ressources, notamment des exemples et l'API.

In [1]:
import cvxpy as cp
import numpy as np

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 = np.array(range(-4,6))   # beta is just -4 through 5.

X = np.random.randn(n,p)
Y = np.dot(X,beta_star) + 0.1*np.random.randn(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 CVXPY

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

objective = cp.Minimize(L)

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

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]:
prob = cp.Problem(objective)
  • Étape 4. Le résoudre!
In [6]:
prob.solve()
Out[6]:
0.10370392166611218
  • Étape 5. Analyser le résultat

Après résolution, la variable prob du solver possède (entre autres) les attributs suivants:

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

value : valeur optimale de L

solver_stats.solve_time : le temps mis à trouver une solution

la variable beta contient notamment la valeur optimal

value : valeur optimale (de la variable précisée)

In [7]:
prob.status
Out[7]:
'optimal'
In [8]:
prob.value
Out[8]:
0.10370392166611218
In [9]:
prob.solver_stats.solve_time
Out[9]:
0.000326303
In [10]:
print("\nThe optimal value is", prob.value)
print("The optimal beta is")
print(beta.value)
print("The norm of the residual is ", cp.norm(X*beta - Y, p=2).value)
The optimal value is 0.10370392166611218
The optimal beta is
[-4.01135908 -2.99415813 -2.02166575 -0.95690415 -0.03439186  0.99785978
  1.96812678  3.00440345  3.9468306   4.97922063]
The norm of the residual is  0.3220309327783794

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, CVXPY 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 [11]:
beta = cp.Variable(p)  # p est la taille de la variable, définie ci-dessus
  • Étape 2. Définir la fonction objectif et les contraintes
In [12]:
L = cp.sum_squares(X*beta - Y) # cost est la fonction ci-dessus 

objective = cp.Minimize(L)
In [13]:
contrainte1 = beta>= np.zeros(p)
In [14]:
ones = np.ones((1,p))
contrainte2 = ones*beta == 1
  • Étape 3. Créer le problème à résoudre avec les contraintes
In [23]:
prob = cp.Problem(objective,constraints=[contrainte1,contrainte2])
  • Étape 4. Le résoudre!
In [16]:
prob.solve()
Out[16]:
874.6956008866483
  • Étape 5. Analyser le résultat
In [17]:
prob.status
Out[17]:
'optimal'
In [18]:
betaHat2 = beta.value
betaHat2
Out[18]:
array([-2.06162116e-15, -1.64618429e-15, -9.51012819e-16, -1.60920304e-15,
       -1.20738874e-15, -1.40756977e-15, -1.34693570e-15, -1.34313678e-15,
       -1.25909279e-15,  1.00000000e+00])
In [19]:
np.sum(betaHat2)
Out[19]:
1.0000000000000018

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 [20]:
lambda1 = prob.constraints[0].dual_value
lambda1
Out[20]:
array([246.04950955, 238.62525452, 135.04977932, 174.3783018 ,
       164.51231503, 180.67554211, 142.75247838,  59.04302607,
        82.42694851,   0.        ])
In [21]:
print("value \t dual of constraint 1")
for i in range(p):
    print("{:2.4f}  \t {:4.2f}".format(betaHat2[i],lambda1[i]))
value 	 dual of constraint 1
-0.0000  	 246.05
-0.0000  	 238.63
-0.0000  	 135.05
-0.0000  	 174.38
-0.0000  	 164.51
-0.0000  	 180.68
-0.0000  	 142.75
-0.0000  	 59.04
-0.0000  	 82.43
1.0000  	 0.00
In [22]:
prob.constraints[1].dual_value
Out[22]:
array([145.98401485])