Algorithmes d'optimisation -- L3 MINT et doubles licences 2019/2020 -- Université Paris-Saclay
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
%matplotlib inline
# la commande suivante agrandit les figures
plt.rcParams['figure.figsize'] = [9.,6.]
def verifier_gradient(f,g,x0):
N = len(x0)
gg = np.zeros(N)
for i in range(N):
eps = 1e-4
e = np.zeros(N)
e[i] = eps
gg[i] = (f(x0+e) - f(x0-e))/(2*eps)
print('erreur numerique dans le calcul du gradient: %g (doit etre petit)' % np.linalg.norm(g(x0)-gg))
$\newcommand{\Rsp}{\mathbb{R}}$ $\newcommand{\sca}[2]{\langle #1|#2\rangle}$ $\newcommand{\eps}{\varepsilon}$ $\newcommand{\proj}{\mathrm{proj}}$ Comme dans le troisième exercice de la feuille de TD consacrée au calcul de projections, on appelle simplexe l'ensemble
$$\Delta = \{ x\in \Rsp_+^n \mid \sum_{1\leq i \leq n} x_i = 1\}.$$Dans un premier temps, on va chercher à calculer la projection d'un point $y\in\Rsp^n$ sur $\Delta$. Pour cela, on admettra les deux résultats suivants, démontrés dans le TD:
Ainsi, pour trouver la projection d'un point $y\in\Rsp^n$ sur le simplexe $\Delta$, il suffit de trouver $\kappa\in\Rsp$ tel que $g(\kappa) = 0$ où l'on a posé $$ g(\kappa) = \left(\sum_{1\leq i\leq n} \max(y_i - \kappa, 0)\right) -1$$
Q1) Soit $y = (0.1,1.5,2.1) \in \Rsp^3$. Écrire la fonction g(kappa)
décrite ci-dessus. Trouver $\kappa$ vérifiant $g(\kappa) = 1$ en utilisant la fonction scipy.optimize.root(g,x0=0,method='anderson').x
.
# <completer>
Q2) En s'inspirant du code de la fonction précédente, écrire une fonction proj_simplexe calculant la projection d'un point $y\in\Rsp^n$ sur $\Delta$. Pour vérifier le bon fonctionnement de cette fonction, calculer calculer $p=$proj_simplexe(y)
pour $y=$np.random.randn(n)
puis vérifier que
où $e_i$ est le $i$ème vecteur de la base canonique.
(Question subsidiaire: pourquoi cette inégalité doit-elle être vraie pour $p = \proj_\Delta(y)$ ? Caractérise-t-elle la projection sur $\Delta$?)
def proj_simplexe(y):
# <completer>
# validation: calcul de produits scalaires
# <completer>
Dans cette partie, il s'agit de résoudre un problème d'optimisation sous contraintes de la forme suivante:
$$ \min_{x\in\Delta} f(x) \hbox{ où } f(x) = \frac{1}{2}\sca{x}{Qx} + \frac{1}{2\eta}(\sca{r}{x}-r_0)^2, $$où $Q \in \mathcal{M}_{n,n}(\Rsp)$ est une matrice symétrique définie positive, $r\in\Rsp$ est un vecteur et $r_0\in\Rsp$ sont donnés. Dans la suite, on fixera $\eta = 10^{-2}$.
Motivation: Ce problème modélise une situation où un investisseur cherche placer un portefeuille en garantissant un certain niveau de rendement $r_0\in\Rsp$ tout en minimisant le risque. Plus précisément, dans ce problème $n$ est le nombre d'actifs (actions, etc.), et la variable $x\in\Rsp^n$ décrit la stratégie d'investissement: $x_i$ décrit la fraction du portefeuille que l'on investit dans l'actif $i$. Ainsi, la contrainte $x\in\Delta$ modélise:
La matrice $Q$ modélise la covariance entre les actifs:
Le problème d'optimisation qu'on regarde consiste donc à trouver une stratégie d'investissement ($x\in \Delta$) cherchant à viser un rendement donné (terme $\frac{1}{2\eta}(\sca{r}{x}-r_0)^2$) tout en minimisant le risque (terme $\sca{x}{Qx}$).
En pratique, on considèrera les données suivantes:
# Le code suivant permet de construire la matrice $Q$ à partir de données réelles.
# Ne pas hésiter à la décommenter et à tester avec d'autres choix d'actifs. Pour cela,
# il faut installer le package yfinance en lançant la ligne suivant
# !pip install yfinance
if False: # mettre True pour changer les dates, les actifs, etc
import yfinance as yf
import numpy as np
stocks = ['NFLX','MSFT','BA','AIR']
closes = np.array([np.array(yf.download(s,'2019-01-01','2020-01-01').Close) for s in stocks]).T
r = (closes[-1,:] - closes[1,:])/closes[1,:] # rendement sur l'année 2019
Q = np.cov(closes - np.mean(closes,0),rowvar=False) # estimation de la covariance entre les valeurs
Q = np.array([[1199.6242199, -225.74269344, 270.42617708, -112.31853678],
[-225.74269344, 224.42514399, -157.75776414, 46.31290714],
[ 270.42617708, -157.75776414, 600.37115079, -28.1665365 ],
[-112.31853678, 46.31290714, -28.1665365, 21.77422792]])
r = np.array([0.20888442, 0.55953316, 0.00602209, 0.2042723])
r0 = 0.7*np.max(r)
eta = 1e-4
Q1) Montrer que la fonction $f$ est strictement convexe et que $\nabla f(x) = Qx + \frac{1}{\eta} (\sca{r}{x}-r_0) r.$
Q2) Écrire une fonction f
et gradf
, et utiliser la fonction verifier_gradient
pour valider l'implémentation.
# <completer>
# vérification du calcul du gradient
verifier_gradient(f,gradf,np.random.rand(len(r)))
Q3) Implémenter l'algorithme du gradient projeté pour résoudre le problème avec les données ci-dessus, i.e.
$$ \begin{cases} x^{(0)} = 0_{\Rsp^n}\\ x^{(k+1)} = \mathrm{proj\_simplexe}(x - \tau \nabla f(x^{(k)})) \end{cases}$$Tracer sur deux figures distinctes:
Pour choisir le paramètre $\tau$ de l'algorithme du gradient projeté, on pourra procéder par tâtonnement: par exemple, choisir $\tau$ de sorte à ce que $k\mapsto f(x^{(k)})$ soit décroissante et que la suite $\|{x^{k} - x^{(k-1)}}\|$ semble tendre vers zéro.
# <completer>
Q1 Montrer que le simplexe $\Delta = \{ x\in\Rsp_+^n\mid \sum_i x_i = 1\}$ peut-être mis sous la forme
$$ \Delta = \{ x\in\Rsp^n \mid \forall 1\leq i\leq \ell, c_i(x) \leq 0 \} $$avec $\ell = n+2$ et
$$ c_i(x) = \begin{cases} -x_i & \hbox{ si } 1\leq i\leq n \\ x_1+\dots+x_n - 1 & \hbox{ si } i=n+1\\ -(x_1+\dots+x_n - 1) & \hbox{ si } i=n+2\end{cases} $$Dans la méthode de pénalisation, le problème d'optimisation sous contraintes
$$ P = \min_{x\in\Delta} f(x) \hbox{ où } f(x) = \frac{1}{2}\sca{x}{Qx} + \frac{1}{2\eta}(\sca{r}{x}-r_0)^2, $$est alors approché par le problème d'optimisation sans contraintes suivant:
$$ P_\eps = \min_{x\in\Rsp^n} f_\eps(x) \hbox{ où } f_\eps(x) = f(x) + \frac{1}{\eps}\sum_{1\leq i\leq \ell} \max(c_i(x),0)^2 $$Q2 Montrer que
$$ f_\eps(x) = f(x) + \frac{1}{\eps}\left(\sum_{1\leq i\leq n} \max(-x_i,0)^2)\right) + \frac{1}{\eps}(x_1+\dots+x_n - 1)^2 $$$$(\nabla f_\eps(x))_i = (\nabla f(x))_i - \frac{2}{\eps} \max(-x_i,0) + \frac{2}{\eps} (x_1+\dots+x_n - 1)$$Coder deux fonctions feps(x)
et gradfeps(x)
, et vérifier leur bon fonctionnement en utilisant verifier_gradient
.
eps = 1e-3
# <completer>
verifier_gradient(feps,gradfeps,10*(np.random.rand(len(r))-.5))
Q3 En utilisant la fonction gradient_armijo
ci-dessous, vérifier le fonctionnement de cette approche pour des valeurs de eps
modérées ($\eps = 10^{-2}$ où $10^{-3}$). Commenter les résultats (respect des contraintes, vitesse de convergence).
def gradient_armijo(f,gradf,x0,err=1e-5,maxiter=20000):
x = x0.copy()
fiter = []
giter = []
k = 0 # nombre d'itérations
while(True):
k = k+1
if k > maxiter: # maximum de 10^6 itérations
print('erreur: nombre maximum d\'itérations atteint')
break
d = -gradf(x)
fiter.append(f(x))
giter.append(np.linalg.norm(d))
if np.linalg.norm(d) <= err:
break
t = 1
m = -np.dot(d,d)
while f(x+t*d) > f(x) + 0.3*t*m:
t = 0.5*t
if k%100==0: # on affiche des informations toute les 100 itérations
print('iteration %d: f=%f, |g|=%f, step=%f' % (k, f(x), np.linalg.norm(d),t))
x = x + t*d
return x,np.array(fiter),np.array(giter)
# <completer>
x
array([0.12850422, 0.48099323, 0.04045281, 0.36423465])
np.sum(x)
1.0141849153713878