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.
import numpy as np
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$. Polynomia
du module numpy.polynomial
.from numpy.polynomial import Polynomial
def valeurPolynome(coefs, x):
p = Polynomial(coefs)
return p(x)
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.coefs = np.array([1, 2, 3])
x1 = 2
print(valeurPolynome(coefs, x1))
x2 = np.array([[1, 2], [3, 4]])
print(valeurPolynome(coefs, x2))
derivPoly(coeff)
qui retourne le tableau de coefficients représentant la dérivée du polynôme défini par coeff
. Testez.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))
%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)
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 ?%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')
nu
au fur et à mesure des itérations ? Quel est l'intérêt de procéder ainsi?%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)
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()
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
.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))
pdep = (1, -5)
puis (-1, -3)
, (-1, -5)
et (2, -4)
.animationDescenteSurface(x,y,z,(1,-5))
animationDescenteSurface(x,y,z,(-1,-3))
animationDescenteSurface(x,y,z,(-1,-5))
animationDescenteSurface(x,y,z,(2,-4))
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 :
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')
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
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')