TME 3 - Descente de gradient

Il s'agit dans ce TME de coder la descente de gradient pour une fonction générique. Nous utiliserons pour cela une classe abstraite OptimFunc, qui contiendra la définition de la fonction à calculer (f), le gradient de cette fonction (grad_f) et la fonction x_random pour générer aléatoirement un point initial.



In [28]:

    
import numpy as np
from numpy import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# fonction pour transformer un vecteur en array si ce n'est pas un array
#, a mettre en première ligne de l'appel de fonctions, comme dans le squelette OptimFunc
def to_array(x):
    if len(x.shape)==1:
        x=x.reshape(1,x.shape[0])
    return x

#Classe abstraite d'une fonction à optimiser. self.dim contient la dimension attendue de l'entrée (pour la génération aléatoire d'un point).
class OptimFunc:
    def __init__(self,dim=2):
        self.dim=dim
    def f(self,x): #fonction a calculer
        x=to_array(x)
        pass
    def grad_f(self,x): #gradient de la fonction
        x=to_array(x)
        pass
    def x_random(self,low=2,high=5):
        return random.random(self.dim)*(high-low)+low

Voici le squelette de la classe pour l'optimisation par descente de gradient. Les paramètres sont les suivants:

eps : epsilon, le pas du gradient; la fonction get_eps() permet de faire un pas adaptatif
optiim_f : objet de type OptimFunc, la fonction à optimiser
max_iter : le nombre d'itération maximum, None s'il ne faut pas le prendre en compte
delta : la différence minimale à avoir entre deux pas pour continuer à itérer

En plus de ces paramètres, les variables suivantes sont utilisées :

i : le numéro de l'itération courante
x : le point optimal courant
log_x, log_f, log_grad : l'historique des points parcourus, de la valeur de la fonction en ces points, de la valeur du gradient en ces points.

Compléter ce qui manque. Utiliser en particuliler np.vstack((log_x,x)) pour empiler au fur et à mesure les résultats dans l'historique.



In [29]:

    
class GradientDescent:
    def __init__(self,optim_f,eps=1e-4,max_iter=5000,delta=1e-6):
        self.eps=eps
        self.optim_f=optim_f
        self.max_iter=max_iter
        self.delta=delta
    def reset(self):
        self.i=0
        self.x = self.optim_f.x_random()
        self.log_x=np.array(self.x)
        self.log_f=np.array(self.optim_f.f(self.x))
        self.log_grad=np.array(self.optim_f.grad_f(self.x))
    def optimize(self,reset=True):
        if reset:
            self.reset()
        while not self.stop():
            self.x = self.x - self.get_eps()*self.optim_f.grad_f(self.x)
            self.log_x=np.vstack((self.log_x,self.x))
            self.log_f=np.vstack((self.log_f,self.optim_f.f(self.x)))
            self.log_grad=np.vstack((self.log_grad,self.optim_f.grad_f(self.x)))
            self.i+=1
    def stop(self):
        return (self.i>2) and (self.max_iter and (self.i>self.max_iter) or (self.delta and np.abs(self.log_f[-1]-self.log_f[-2]))<self.delta)
    def get_eps(self):
        return self.eps

        
        
class Rosenbrock(OptimFunc):   
    def f(self,x): #fonction a calculer
        x=to_array(x)
        return 100*(x[:,1]-x[:,0]**2)**2 + (1-x[:,0])**2
    def grad_f(self,x): #gradient de la fonction
        x=to_array(x)
        x1=x[:,0]
        x2=x[:,1]
        d1=200*(x2-x1**2)*(-2*x1)-2*(2-2*x1)
        d2=200*(x2-x1**2)
        return np.array([d1,d2]).T
    
class FF1(OptimFunc):   
    def f(self,x): #fonction a calculer
        x=to_array(x)
        return x*np.cos(x)
    def grad_f(self,x): #gradient de la fonction
        x=to_array(x)
        return np.cos(x)-x*(np.sin(x))   

class FF2(OptimFunc):   
    def f(self,x): #fonction a calculer
        l=[]
        for i in x:
            if i==0:
                l.append(0)
            else:
                l.append(-1*np.log(i) + (i)**2)
        return np.asarray(l)
    
    def grad_f(self,x): #gradient de la fonction
        x=to_array(x)
        return -1/x + 2*x

Mise en pratique

Tester votre implémentation sur 3 fonctions :

en 1d sur la fonction f(x)=x.cos(x)
en 1d sur la fonction -log(x)+x^2
en 2d sur la fonction Rosenbrock (ou banana), définie par : f(x_1,x_2)= 100*(x_2-x_1^2)^2 + (1-x_1)^2 .

Tracer les courbes d'évolution de f(x) en fonction du nombre d'itérations. Tracer (en 3d ou en 2d) la surface de la fonction, et la trajectoire de l'optimisation. Vous pouvez utilisez le code suivant pour la visualisation 3d.

Descente de gradient sur la fonction Rosenbrock



In [5]:

    
#Crétion de la fonction Rosenbrock
rosen = Rosenbrock()

#Quadrillage pour evaluer la fonction f rosenbrock.
m=5
x=np.arange(-m,m,2*m/20.)
y=np.arange(-m,m,2*m/20.)
xx,yy=np.meshgrid(x,y)
grid=np.c_[xx.ravel(),yy.ravel()]
z=rosen.f(grid).reshape(xx.shape)

#Création du gradient et optimisation de celui ci
grd_rosen = GradientDescent(rosen)
grd_rosen.optimize()

#Création de la figure
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, yy, z, rstride=1, cstride=1, cmap=cm.gist_rainbow, linewidth=0, antialiased=False)
fig.colorbar(surf)


ax.plot(grd_rosen.log_x[:,0], grd_rosen.log_x[:,1], grd_rosen.log_f.ravel(),color='black')
titre="Gradient_descent_Rosenbrock"
plt.title(titre)
savefig(titre)
plt.show()

Descente de gradient sur la fonction f(x)=x.cos(x)



In [8]:

    
#création de la fonction
ft = FF1(1)

#création de x et y
x = np.arange(-6,6,0.1)
f_x = ft.f(x).ravel()


#création du gradient et optimisation
grd = GradientDescent(ft,eps=1e-4,max_iter=200000,delta=1e-7)  
grd.optimize()

plt.plot(x,f_x)
plt.plot(grd.log_x, grd.log_f, color='red')
titre="Gradient_descent_xcos(x)"
plt.title(titre)
savefig(titre)
plt.show()

Cheminements des gradients

On remarque que le gradient se dirige vers un unique minimum local ici il y en a 3.



In [6]:

    
plt.plot(x,f_x)
for i in range(30):
    grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
    grd.optimize()
    plt.plot(grd.log_x, grd.log_f, color='red')
    
titre="Cheminement_gradient"
plt.title(titre)
savefig(titre)
plt.show()

Ordonnée du gradient en fonction du nombre d'itérations.

On remarque donc que la fonction possède 3 minimums locaux en y=-0.5 y=-3.5 et y=-6



In [22]:

    
for i in range(30):
    grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
    grd.optimize()
    plt.plot(ft.f(grd.log_x))
    
titre="Gradient_descent_ordonnee_xlog(x)"
plt.title(titre)
savefig(titre)
plt.show()

Abcisse du gradient en fonction du nombre d'itérations.

On remarque donc que la fonction possède 3 minimums locaux en x=-6 x=-1 et x=3.2



In [17]:

    
for i in range(30):
    grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
    grd.optimize()
    plt.plot(grd.log_x)
titre="Gradient_descent_abscisse_xlog(x)"
plt.title(titre)
savefig(titre)
plt.show()

Descente de gradient sur la fonction -log(x)+x^2



In [30]:

    
#création de la fonction
ft = FF2(1)

#création de x et y
x = np.arange(0.2,5,0.1)
f_x = ft.f(x).ravel()


#création du gradient et optimisation
grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
grd.optimize()

plt.plot(x,f_x)
plt.plot(grd.log_x, grd.log_f, color='red')
titre="Gradient_descent_-log(x)+x^2"
plt.title(titre)
savefig(titre)

plt.show()

Ordonnée du gradient en fonction du nombre d'itérations.

On remarque donc que la fonction possède 1 minimums locaux en y=1



In [31]:

    
for i in range(30):
    grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
    grd.optimize()
    plt.plot(ft.f(grd.log_x))
titre="Gradient_descent_ordonnee_-log(x)+x^2"
plt.title(titre)
savefig(titre)
plt.show()

Abcisse du gradient en fonction du nombre d'itérations.

On remarque donc que la fonction possède 1 minimums local en x=1



In [32]:

    
for i in range(30):
    grd = GradientDescent(ft,eps=1e-4,max_iter=20000,delta=1e-7)  
    grd.optimize()
    plt.plot(grd.log_x)
titre="Gradient_descent_abscisse_-log(x)+x^2"
plt.title(titre)
savefig(titre)
plt.show()



In [ ]: