115 KiB
115 KiB
In [1]:
import numpy as np
import matplotlib.pyplot as plt
Question 1¶
In [ ]:
def gradientquad(Q,c,p,x0,maxiter=100,choix='exact'):
fobj = np.zeros(maxiter)
x = x0.copy()
valp, vecp = np.linalg.eig(Q)
L = np.max(valp)
# Si le choix est "constant", alpha ne dépend pas de l'itéré, on le calcule une seule fois.
alpha = 1/L
for k in range(maxiter):
# Gradient de fonction quadratique
d = -(Q@x-c)
if choix == "exact":
# Dépend de l'itéré, on le calcule à chaque itéré
alpha = (d@d)/(d@Q@d)
x = x + alpha*d
fobj[k] = 0.5*x@Q@x - c@x + p
if np.linalg.norm(d) <= 1e-6:
break
return x, fobj, k
Expérience 1¶
In [3]:
Q = np.array([[10.0,-6.0],[-6.0,10.0]])
c = np.array([4.0,-3.0])
p = 0
x0 = np.array([1.0,1.0])
# Méthode exacte
xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=25,choix='exact')
plt.plot(fobj, 'b', label="Gradient avec pas exact")
# Méthode constante
xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=25,choix='constant')
plt.plot(fobj, 'r', label="Gradient avec pas constant")
plt.legend()
Out[3]:
Expérience 2¶
In [4]:
from scipy.linalg import hilbert
In [5]:
nbiter_exact = []
nbiter_constant = []
for n in range(3,101):
Q = hilbert(n)
c = np.zeros(n)
p = 0
x0 = np.ones(n)
# Méthode exacte
xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=10**6,choix='exact')
nbiter_exact.append(nbiter)
# Méthode constante
xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=10**6,choix='constant')
nbiter_constant.append(nbiter)
plt.plot(nbiter_exact, label="Pas exact")
plt.plot(nbiter_constant, label="Pas constant")
plt.xlabel("Dimension n")
plt.ylabel("Nombre d'itérations")
Out[5]:
Question 2¶
Implémenter les fonctions¶
In [10]:
def f_grad_conv(x, c=10.0):
f = 2*np.log(np.cosh((c/2)*(x[0]+2*x[1]))) + 2*np.log(np.cosh((c/2)*(2*x[0]-x[1])))
grad = np.array([
c*np.tanh((c/2)*(x[0]+2*x[1])) + 2*c*np.tanh((c/2)*(2*x[0]-x[1])),
2*c*np.tanh((c/2)*(x[0]+2*x[1])) - c*np.tanh((c/2)*(2*x[0]-x[1]))
])
L = 2.5*c**2
return f, grad, L
In [11]:
def f_grad_strong(x, c=10.0, mu=10**2):
f, grad, L = f_grad_conv(x, c)
f += (mu/2)*(x[0]**2 + x[1]**2)
grad += mu*x
L += mu
return f, grad, L
Implémenter la fonction methgradient¶
In [43]:
def methgradient(fct, x0, maxiter=100, tol=1e-9):
x = x0.copy()
fobj = np.zeros(maxiter)
for k in range(maxiter):
f, grad, L = fct(x)
d = grad
alpha = 1/L
x = x - alpha*d
f, grad, L = fct(x)
fobj[k] = f
if np.linalg.norm(grad) <= tol:
break
return x, fobj[:k]
Expérience¶
In [44]:
x0 = np.array([1.0,1.0])
xopt, fobj = methgradient(f_grad_conv, x0, maxiter=70)
plt.semilogy(fobj, 'b', label="Fonction convexe")
xopt, fobj = methgradient(f_grad_strong, x0, maxiter=70)
plt.semilogy(fobj, 'r', label="Fonction fortement convexe")
plt.xticks(np.arange(0,18,2))
plt.legend(loc='upper right')
Out[44]:
Question 3¶
Implémenter la fonction¶
Note: La fonction n'est pas L-smooth parce qu'on peut pas trouver une fonction quadratique avec un seul L tel que cette fonction est au dessus de $x^4$. La courbature de $x^4$ en $\infty$ est aussi $\infty$.
In [45]:
def fct(x):
f = (x[0]+x[1])**4 - 2*(x[0]+x[1])**2 + 1
g1 = 4*(x[0]+x[1])**3 - 4*(x[0]+x[1])
g = np.array([g1, g1])
return f, g
Implémenter les deux conditions de Wolfe¶
In [104]:
def w1(fct, x, d, alpha, beta1):
f0, g0 = fct(x)
f1, _ = fct(x+alpha*d)
return f1 <= f0 + beta1*alpha*(g0@d)
def w2(fct, x, d, alpha, beta2):
_, g0 = fct(x)
_, g1 = fct(x+alpha*d)
return d@g1 >= beta2*(d@g0)
Implémenter la bissection de Wolfe¶
In [105]:
def wolfebissection(fct, x, d, alpha0=1.0, beta1=1e-4, beta2=0.9):
aleft = 0
aright = np.inf
alpha = alpha0
while True:
if w1(fct, x, d, alpha, beta1) and w2(fct, x, d, alpha, beta2):
return alpha
if not w1(fct, x, d, alpha, beta1):
# On bouge la borne de droite vers la gauche
aright = alpha
alpha = (aleft + aright)/2
elif not w2(fct, x, d, alpha, beta2):
# On bouge la borne de gauche vers la droite
aleft = alpha
if aright < np.inf:
alpha = (aleft + aright)/2
else:
alpha *= 2
Implémenter la méthode du gradient avec recherche de Wolfe¶
In [106]:
def methgradient(fct, x0, maxiter=100, tol=1e-6, alpha0=1e-6):
x = x0.copy()
fobj = np.zeros(maxiter)
alpha = alpha0
for k in range(maxiter):
f, g = fct(x)
d = -g
alpha = wolfebissection(fct, x, d, alpha)
x = x + alpha*d
fobj[k] = f
if np.linalg.norm(g) <= tol:
break
return x, fobj[:k]
Exemple numérique¶
In [107]:
x0 = np.array([1.0,1.0])
maxiter = 30
# Importance de alpha0!
xopt, fobj = methgradient(fct, x0, maxiter=maxiter, alpha0=0.05)
plt.plot(fobj, 'b', label="Gradient avec Wolfe")
plt.legend(loc="upper right")
print(xopt)
Bonus: Méthode gradient accéléré¶
In [108]:
def methgradientaccel(fct, x0, maxiter=100, tol=1e-6, alpha0=1e-6):
x = x0.copy()
fobj = np.zeros(maxiter)
xprev = x
alpha = alpha0
for k in range(maxiter):
beta = k/(k+3) # Paul Tseng (version plus simple)
y = x + beta*(x - xprev)
f, g = fct(y)
d = -g
alpha = wolfebissection(fct, y, d, alpha)
xprev = x.copy()
x = y + alpha*d
f, g = fct(x)
fobj[k] = f
if np.linalg.norm(g) <= tol:
break
return x, fobj[:k]
In [ ]:
x0 = np.array([1.0,1.0])
maxiter = 30
# Importance de alpha0!
xopt, fobj = methgradientaccel(fct, x0, maxiter=maxiter, alpha0=0.05)
plt.plot(fobj, 'b', label="Gradient avec Wolfe")
plt.legend(loc="upper right")
print(xopt)