Université de Bordeaux, Master Mention Informatique

Analyse, classification et indexation des données: feuille 2

Descente de gradient

Avant de commencer :

La méthode de la descente de gradient permet de trouver un minimum local d’une fonction de ${\mathbb R}^n$ dans ${\mathbb R}$ . L’objectif de ce TD est de comprendre cette méthode pour pouvoir l’utiliser par la suite dans différentes applications où on cherche à minimiser un critère.

Comme précédemment nous utilisons entre autres numpy pour manipuler vecteurs et matrices.

In [1]:
import numpy as np

Exercice 1. Polynôme : calcul et dérivée

  1. Ecrivez une fonction python valeurPolynome (coeff, x) qui calcule la valeur du polynôme défini par le vecteur de coefficients coeff pour la valeur x. Les coefficients sont rangés dans l'ordre des puissances croissantes. Par exemple le vecteur $[1, 2, 3]$ définit le polynôme $3x^2+2x+1$.
    Utilisez pour cela la classe Polynomia du module numpy.polynomial.
In [2]:
from numpy.polynomial import Polynomial
def valeurPolynome(coefs, x):
    p = Polynomial(coefs)
    return p(x)
  1. Testez cette fonction pour différentes valeurs de x et pour un polynôme dont vous choisirez les coefficients. Remarquez que le paramètre x peut lui même être un vecteur ou une matrice.
In [3]:
coefs = np.array([1, 2, 3])
x1 = 2
print(valeurPolynome(coefs, x1))
x2 = np.array([[1, 2], [3, 4]])
print(valeurPolynome(coefs, x2))
17.0
[[ 6. 17.]
 [34. 57.]]
  1. La dérivée d’un polynôme est elle même un polynôme. Ecrivez une fonction python derivPoly(coeff) qui retourne le tableau de coefficients représentant la dérivée du polynôme défini par coeff. Testez.
In [4]:
def derivPoly(coefs):
    p = Polynomial(coefs)
    return p.deriv().coef

#On teste le résultat :
coefs = np.array([1, 2, 3])
p = Polynomial(coefs)
print(p.coef)
print(derivPoly(coefs))
[1. 2. 3.]
[2. 6.]

Exercice 2. Descente de gradient pour trouver le minimum d’une parabole

  1. Exécutez le script python suivant :
In [5]:
%matplotlib 
import matplotlib.pyplot as plt

def parabole(x):
    return x * x

x = np.linspace(-300,300,600)
y = [parabole(i) for i in x]

epsilon = 0.1
nu = 0.2

plt.plot(x,y,'g')

x0=250
xcurrent = x0 - 2*x0
xprec = x0

plt.plot(x0,parabole(x0),'ok')


while(abs(xprec - xcurrent) > epsilon):
    xprec  = xcurrent
    xcurrent = xprec - nu*2*xprec
    plt.plot(xcurrent,parabole(xcurrent),'ob')
    plt.pause(1)
    
plt.plot(xcurrent,parabole(xcurrent),'xr')
print(xcurrent)
Using matplotlib backend: TkAgg
-0.11754624614399992
  1. Retrouvez dans le code les étapes de l’algorithme de descente du gradient.
  2. Testez en modifiant l’initialisation de la position du minimum.
  3. Comment peut-on obtenir une estimation plus précise de la position du minimum?
  4. Testez en modifiant la variable nu par une autre valeur constante inférieure à 1, par exemple nu = 0.6. Que se passe-t-il si la valeur de nu est fixée à 1 ? Est-ce spécifique à la fonction choisie ?
  5. Comparez la version suivante du script avec la précédente
In [6]:
%matplotlib
import matplotlib.pyplot as plt

def parabole(x):
    return x * x

x = np.linspace(-300,300,600)
y = [parabole(i) for i in x]

epsilon = 0.1
attenuation = 3

plt.plot(x,y,'g')

x0=250
xcurrent = x0 - 2*x0
xprec = x0

plt.plot(x0,parabole(x0),'ok')

while(abs(xprec - xcurrent) > epsilon):
    xprec  = xcurrent
    nu = 1/np.log(attenuation)
    attenuation = attenuation+1
    xcurrent = xprec - nu*2*xprec
    plt.plot(xcurrent,parabole(xcurrent),'ob')
    plt.pause(1)
plt.plot(xcurrent,parabole(xcurrent),'xr')
Using matplotlib backend: TkAgg
Out[6]:
[<matplotlib.lines.Line2D at 0x7f2869592ac0>]
  1. Comment évolue la variable nu au fur et à mesure des itérations ? Quel est l'intérêt de procéder ainsi?

Exercice 3. Descente de gradient pour trouver le minimum d’un polynôme

  1. On choisit comme exemple le polynôme $30 − 61x + 41x^2 − 11x^3 + x^4$ . Ecrivez un script python qui affiche ce polynôme entre 0 et 6 et qui montre les étapes de la descente de gradient à partir de la valeur initiale 5.
In [7]:
%matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0,6,60)
nu = .01
epsilon = 0.01 
coefpoly = np.array([30, -61, 41, -11, 1]) 

y = [valeurPolynome(coefpoly, xi) for xi in x] 

derive = derivPoly(coefpoly)

plt.plot(x,y)

x0 = 5
plt.plot(x0, valeurPolynome(coefpoly, x0), 'xk')
xprec = x0
xcurrent = xprec - nu*valeurPolynome(derive, xprec)

while(np.abs(xprec - xcurrent) > epsilon):
    plt.plot(xcurrent, valeurPolynome(coefpoly, xcurrent), 'xr')
    xprec  = xcurrent
    xcurrent = xprec - nu*valeurPolynome(derive, xprec)
    plt.pause(0.1)
    
Using matplotlib backend: TkAgg
  1. Changez la valeur initiale à 0.
  2. Faites varier les paramètres de la descente de gradient (atténuation, test d’arrêt) et commentez le comportement de l’algorithme.

Exercice 4. Descente de gradient pour trouver le minimum d’une surface

  1. On choisit comme exemple la fonction définie par $f(x,y) = (x−1)(x−2)+(y+3)(y+4)$. Visualisez la surface correspondante avec $x$ et $y$ variant entre $−8$ et $8$ en complétant le code ci-dessous. Quelles sont les coordonnées du minimum ?
In [8]:
from mpl_toolkits.mplot3d import Axes3D

def f(x,y):
    return (x - 1) * (x - 2) + (y + 3) * (y + 4)

x = np.linspace(-8, 8, 160)
y = np.linspace(-8, 8, 160)
x, y = np.meshgrid(x, y)

z = f(x, y)

ax = Axes3D(plt.figure())

ax.plot_surface(x, y, z) 
plt.show()
  1. Ecrivez une fonction python animationDescenteSurface(pdep) qui affiche cette surface et qui montre les étapes de la descente de gradient à partir du point initial pdep donné comme vecteur de deux coordonnées. Pour afficher un point 3D vous utiliserez la fonction scatter de Axes3D.
In [9]:
from numpy.linalg import norm

def gradientSurfaceTest(x, y):
    return 2*x - 3, 2*y + 7

def animationDescenteSurface(pdep):
   
    ax = Axes3D(plt.figure())
    
    # affichage point initial
    ax.scatter(pdep[0],pdep[1], f(pdep[0], pdep[1]))

    #affichage surface
#    x = np.linspace(-8, 8, 160)
#    y = np.linspace(-8, 8, 160)
#    x, y = np.meshgrid(x, y)
#    z = f(x, y)
#    ax.plot_surface(x, y, z) 

    xprec = pdep[0]
    yprec = pdep[1]

    epsilon = 1 #0.1

    nu = 0.1 #0.05

    dx,dy = gradientSurfaceTest(xprec,yprec)

    dist = norm([dx,dy]) 
    xcurrent = xprec - nu*dx
    ycurrent = yprec - nu*dy

    while(dist > epsilon):
        ax.scatter(xcurrent, ycurrent, f(xcurrent, ycurrent))
        xprec = xcurrent
        yprec = ycurrent
    
        dx,dy = gradientSurfaceTest(xprec,yprec)
    
        dist = norm([dx,dy]) 
    
        xcurrent = xprec-nu*dx
        ycurrent = yprec-nu*dy
    
        plt.pause(1);
    ax.scatter(xcurrent, ycurrent, f(xcurrent, ycurrent), c='black')
    print(xcurrent)
    print(ycurrent)

animationDescenteSurface((1,2))
1.465640261632
-3.122042877952
  1. Testez avec comme point initial pdep = (1, -5) puis (-1, -3), (-1, -5) et (2, -4).
In [10]:
animationDescenteSurface(x,y,z,(1,-5))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_15487/1444657849.py in <module>
----> 1 animationDescenteSurface(x,y,z,(1,-5))

TypeError: animationDescenteSurface() takes 1 positional argument but 4 were given
In [ ]:
animationDescenteSurface(x,y,z,(-1,-3))
In [ ]:
animationDescenteSurface(x,y,z,(-1,-5))
In [ ]:
animationDescenteSurface(x,y,z,(2,-4))

Exercice 5. Application : régression linéaire

On veut calculer la droite de régression linéaire correspondant à un nuage de points 2D $\{(x_1, y_1), \cdots , (x_m, y_m)\}$, c'est-à-dire la droite $y = \theta_1x + \theta_0$ qui minimise le critère suivant : $$ \frac 1 m \sum \left( \theta_0 + \theta_1 x_i - y_i\right)^2. $$

Autrement dit, on cherche à minimiser la moyenne des distances verticales entre chaque point du nuage et la droite.
Utilisez une descente de gradient (en 2D) pour trouver les caractéristiques $\theta_0$ et $\theta_1$ de la droite de régression linéaire correspondant à un nuage de points. Pour générer un nuage de point, vous pouvez utiliser le code suivant :

In [ ]:
from random import random
m = 40
sizeNoise = 10
x = [random()*50 + 5 for i in range(m)]
noise = [random() * sizeNoise for i in range(m)] # ou np.rand(m,1) * sizeNoise
pente = 0.8
b =  20
y = [b + pente*x[i] + noise[i] for i in range(m)] 

plt.plot(x,y,'bx')
In [ ]:
def gradientJ (theta, x, y):
    theta0 = theta[0]
    theta1 = theta[1]
    m = len(x)
    
    t0 = [theta0 + theta1 * x[i] - y[i] for i in range(m)]
    t1 = [x[i] * (theta0 + theta1 * x[i] - y[i]) for i in range(m)]

    dtheta0 = np.sum(t0)#theta0 + theta1 * x - y)
    dtheta1 = np.sum(t1)#x * (theta0 + theta1*x - y))
    
    dtheta0 = 2*dtheta0 / m
    dtheta1 = 2*dtheta1 / m
    return dtheta0, dtheta1
In [ ]:
from random import random

m = 40
sizeNoise = 10
x = [random()*50 + 5 for i in range(m)]
noise = [random() * sizeNoise for i in range(m)]#np.rand(m,1) * sizeNoise
pente = 0.8
b =  20
y = [b + pente*x[i] + noise[i] for i in range(m)] 


# descente de gradient

theta = (10, 1)
epsilon = 0.1
nu = 0.0005

dtheta0, dtheta1 = gradientJ (theta, x, y)

dist = norm([dtheta0,dtheta1]) 
  
while(dist > epsilon):
   
    theta = theta[0] - nu * dtheta0, theta[1] - nu * dtheta1
    
    dtheta0, dtheta1 = gradientJ (theta, x, y)
    
    dist = norm([dtheta0,dtheta1])

print(theta)


plt.plot(x, y, 'o')
xd = [0, 60]
yd = [theta[0], 60*theta[1]+ theta[0]]
plt.plot(xd , yd, 'r')
In [ ]: