import numpy as np
import matplotlib.pyplot as plt
# comme en TP, on fournit des fonctions permettant de vérifier
# les calculs de gradient et/ou de hessienne
def verifier_derivee(f,gradf,x0):
N = len(x0)
gg = np.zeros(N)
for i in range(N):
eps = 1e-5
e = np.zeros(N)
e[i] = eps
gg[i] = (f(x0+e) - f(x0-e))/(2*eps)
print('erreur numérique dans le calcul du gradient: %g (doit être petit)' % np.linalg.norm(gradf(x0)-gg))
# effectue la descente de gradient avec pas d'Armijo:
# Arguments: f = fonction à minimiser
# gradf = fonction calculant le gradient de d
# x0 = point de départ
# Retourne: x, F où
# x est la solution trouvée
# F = [f(x^(0)), f(x^(1)), ...]
def gradient_armijo(f,gradf,x0):
niter = 200
x = x0.copy()
F = np.zeros(niter)
for k in range(niter):
d = -gradf(x)
m = -np.dot(d,d)
t = 1
while f(x+t*d) > f(x) + 0.3*t*m:
t = 0.5*t
F[k] = f(x)
x = x + t*d
return x,F
Dans cette partie, on considère la fonction $g$ donnée par la formule (3), i.e.
$$ g(t) = a t^2 + b t + \sqrt{\varepsilon + (c+t)^2} \hbox{ où } \varepsilon>0, a\geq 0 $$QN1 Écrire trois fonctions g
, gprime
et gseconde
prenant en argument a,b,c,eps,t
et retournant respectivement $g(t)$, $g'(t)$ et $g''(t)$.
def g(a,b,c,eps,t):
# à compléter
# définir gprime et gseconde
## vérification du calcul de gprime et gseconde
a=0.1; b=2; c=-3; eps=0.01
verifier_derivee(lambda t: g(a,b,c,eps,t),lambda t: gprime(a,b,c,eps,t),np.random.rand(1))
verifier_derivee(lambda t: gprime(a,b,c,eps,t),lambda t: gseconde(a,b,c,eps,t),np.random.rand(1))
QN2 Écrire une fonction newton_pur(a,b,c,eps)
qui implémente la méthode de Newton pure donnée par (2):
$$ \begin{cases}
t^{(0)} = 0,\\
t^{(k+1)} = t^{(k)} - g'(t^{(k)})/g''(t^{(k)}).
\end{cases}$$
et retourne la valeur après cinq itérations, i.e. $t^{(5)}$. On utilisera les fonctions gprime(a,b,c,eps,t)
pour calculer $g'(t)$.
def newton_pur(a, b, c, eps):
# à compléter
## vérification du fonctionnement de la fonction Newton pur
print("cas 1: solution trouvée %g, solution attendue -5.0003905411" % newton_pur(0.1,2,-3,0.01))
print("cas 2: solution trouvée %g, solution attendue -995.0000000024029" % newton_pur(0.1,200,-25,0.001))
Il s'agit de résoudre le problème de minimisation (4)
$$ \min_{x\in \mathbb{R}^d} f(x) \hbox{ où } f(x) = \frac{\mu}{2} \|Ax -
y\|^2 + \sum_{1\leq i\leq d} s(x_i) \hbox{ et } s(t) = \sqrt{\varepsilon + t^2}.$$
où $A,y,\mu$ et $\varepsilon$ sont les suivants. Pour vérifier votre code, le minimiseur
du problème est stocké dans xsol
et le minimum dans fmin
# Données du problème
A = np.array([[1,0,1,0,1,0,0,0],
[1,0,0,1,0,1,0,0],
[0,1,1,0,0,0,1,0],
[0,1,0,1,1,0,0,1]])
y = np.array([-1,1,-4,2])
n,d = A.shape
eps = 1e-4
mu = 100
# solution du problème, pour vérification!
fmin = 6.02030542092432
xsol = [0.0000017522, -0.000002031, -1.5035155249, 0.9999950239,
0.5035136323, 0.0000016119, -2.486482525, 0.4864954868]
QN3 Écrire deux fonctions f
, gradf
prenant en argument x
et retournant respectivement $f(x)$ et $\nabla f(x)$.
def f(x):
# à compléter
# définir gradf
## vérification du gradient de f
verifier_derivee(f,gradf,np.random.rand(d))
QN4 Utiliser la fonction gradient_armijo
fournie en préambule pour minimiser $f$ en partant
de x0=(0,...,0)
. Tracer sur une courbe l'évolution de la valeur de la fonction au cours des itérations.
# à compléter
Dans la suite, on mets en oeuvre un algorithme de descente coordonnée par coordonnée.
QN5 Écrire une fonction $T(x,i)$ calculant le minimum du problème (5) $$ \min_{t\in \mathbb{R}} f(x + t e_i),$$ où $e_i$ est le $i$ème vecteur de la base canonique. Pour ce faire,
newton_pure
pour calculer le minimiseur de $g$.# à compléter
L'algorithme de descente coordonnée par coordonnée s'écrit sous la forme suivante. Pour $k\in \mathbb{N}$, on notera $i_k = 1 + \mathrm{mod}(k,d) \in \{1,\dots,d\}$, où $\mathrm{mod}(k,d)$ est le reste de la division euclidienne de $k$ par $d$. On considère les itérées $x^{(k)}$ définies par
$$ x^{(0)}\in\mathbb{R}^d,\quad x^{(k+1)} = x^{(k)} + T(x^{(k)}, i_k) e_{i_k} $$QN6 Tracer $f(x^{(k)}) - $ fmin
en échelle logarithmique, où $1\leq k\leq 2000$ et la
suite $x^{(k)}$ est définie par la formule ci-dessus. On rappelle que pour calculer le reste de la division euclidienne de $k$ par $d$ on peut écrire k % d
.
# à compléter