TP1 : Méthodes de descente de gradient et algorithmes de Newton¶
Menu¶
- Préambule : définition des fonctions
1.1 Étude la fonction de Rosenbrock
1.2 Étude d'une fonction quadratique - Gradient à pas fixe et à pas normalisé
2.1 Implémentation des méthodes
2.2 Application à la fonction de Rosenbrock
2.3 Application à une fonction quadratique - Règle d'Armijo et règle de Wolfe
3.1 Implémentation des règles et descentes de gradient associées
3.2 Application à la fonction de Rosenbrock
3.3 Application à une fonction quadratique - Gradient conjugué
4.1 Implémentation des méthodes
4.2 Application à la fonction de Rosenbrock
4.3 Application à une fonction quadratique - Problème de Lennard-Jones
5.1 Le potentiel de van der Walls
5.2 Application pour $N=4$ et $N=13$ - Méthode Levenberg-Marquardt, Gauss-Newton et régression non-linéaire
6.1 Préambule : régression non-linéaire
6.2 Implémentation de la méthode de Gauss-Newton et de Levenberg-Marquardt
6.3 Complément : test sur le problème de régression linéaire
Préambule : définition des fonctions¶
On commence par charger les paquets qui seront utiles au TP,
In [1]:
import numpy as np
import numpy.random as rnd
import scipy
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
On définit ensuite les fonctions de Rosenbrock et quadratiques qui seront utilisées dans tout le TP. On en profite pour définir leur gradient et leur hessienne.
In [2]:
def mk_quad(m,M,ndim=2):
def quad(x):
y = np.copy(np.asarray(x))
scal = np.ones(ndim)
scal[0] = m
scal[1] = M
y *= scal
return np.sum(y**2)
def quad_grad(x):
y = np.asarray(x)
scal = np.ones(ndim)
scal[0] = m
scal[1] = M
return 2*scal*y
def quad_hessian(x):
scaling = np.ones(ndim)
scaling[0] = m
scaling[1] = M
return 2*np.diag(scaling)
return quad, quad_grad, quad_hessian
def rosenbrock(x):
y = np.asarray(x)
return np.sum((y[0] - 1)**2 + 100*(y[1] - y[0]**2)**2)
def rosenbrock_grad(x):
y = np.asarray(x)
grad = np.zeros_like(y)
grad[0] = 400*y[0]*(y[0]**2-y[1]) + 2*(y[0]-1)
grad[1] = 200*(y[1]-y[0]**2)
return grad
def rosenbrock_hessian_(x):
y = np.asarray(x)
return np.array((
(1 - 4*100*y[1] + 12*100*y[0]**2, -4*y[0]*100),
( -4*100*y[0], 2*100),
))
Étude de la fonction de Rosenbrock¶
On étudie la fonction de Rosenbrock. Cette fonction est très souvent utilisée pour tester des algorithmes d'optimisation. En effet, une fois la "vallée" trouvée, il est assez compliqué pour un algorithme de suivre cette "vallée" jusqu'au point $(1,1)$.
In [3]:
rosenbrock_display = lambda x,y: (1-x)**2+100*(x**2-y)**2
nx = 1000 # number of discretization points along the x-axis
ny = 1000 # number of discretization points along the y-axis
a=-0.5; b=2. # extreme points in the x-axis
c=-1.5; d=4 # extreme points in the y-axis
X,Y = np.meshgrid(np.linspace(a,b,nx), np.linspace(c,d,ny))
Z = rosenbrock_display(X,Y)
CS = plt.contour(X,Y,Z,np.logspace(-0.5,3.5,20,base=10))
plt.title(r'$\mathrm{Rosenbrock \ function: } f(x,y)=(x-1)^2+100(x^2 - y)^2$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Étude d'une fonction quadratique¶
On présente une étude similaire pour la fonction quadratique en deux dimensions. Pour obtenir un comportement de type "vallée étroite" comme on a pu le voir sur la fonction de Rosenbrock, on doit augmenter le rapport $\frac{M}{m}$.
In [4]:
quad_display = lambda x,y,M,m: M*x**2 + m*y**2
M = 10
m = 1
nx = 1000 # number of discretization points along the x-axis
ny = 1000 # number of discretization points along the y-axis
a=-2; b=2. # extreme points in the x-axis
c=-2; d=2 # extreme points in the y-axis
X,Y = np.meshgrid(np.linspace(a,b,nx), np.linspace(c,d,ny))
Z = quad_display(X,Y,M,m)
CS = plt.contour(X,Y,Z,np.linspace(0,10,15))
plt.title(r'$\mathrm{Quadratic \ function: } f(x,y)=Mx^2 + my^2$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Gradient à pas fixe et gradient à pas normalisé¶
Implémentation des méthodes¶
On définit maintenant la descente de gradient à pas fixe.
In [5]:
def steepest_descent_constant_step(f, grad, x0, iterations, error_point, error_grad, h):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = np.asarray(x0)
x_old = x
grad_x = grad(x)
for i in xrange(iterations):
x = x - h*grad(x)
grad_x = grad(x)
f_x = f(x)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
Application à la fonction de Rosenbrock¶
On fixe les paramètres pour la fonction de Rosenbrock et on applique la descente de gradient à pas fixe.
In [6]:
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
h = 10**-3
iterations = 10000
result = steepest_descent_constant_step(f, grad, x0, iterations, error_point, error_grad, h)
In [7]:
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: fixed \ step \ size \ gradient}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Ce comportement est attendu. En effet, une fois que le gradient à pas fixe se retrouve dans la vallée étroite, le gradient est quasiment nul et la descente vers l'argument minimum de la fonction de Rosenbrock est très lente. On reprend maintenant le code de la descente de gradient à pas fixe et on introduit la descente de gradient à pas normalisé.
In [8]:
def steepest_descent_normalized_step(f, grad, x0, iterations, error_point, error_grad, h):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = np.asarray(x0)
x_old = x
grad_x = grad(x)
for i in xrange(iterations):
x = x - h*grad(x)/np.linalg.norm(grad(x))
grad_x = grad(x)
f_x = f(x)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
In [9]:
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
h = 10**-2
iterations = 10000
result = steepest_descent_normalized_step(f, grad, x0, iterations, error_point, error_grad, h)
In [10]:
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: normalized \ step \ size \ gradient}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
La normalisation permet à l'algorithme de ne pas exploser pour des valeurs élevées du pas, la quantité de déplacement est toujours la même. Si cela peut avoir un avantage, on a par contre plus de mal à obtenir des résultats précis. En effet, puisque la quantité de déplacement est toujours la même on ne peut pas approcher le pas à moins de $10^{-2}$. Cette valeur est atteinte dans notre cas (voir $\texttt{point error}$).
Application à une fonction quadratique¶
On présente maintenant une étude similaire pour la fonctionnelle quadratique introduite en première partie.
In [11]:
M = 5.
m = 1.
n = 3
(f,grad,hess) = mk_quad(M,m,n)
x0 = np.array([5.,5.,5.])
error_point = 10**-10
error_grad = 10**-10
h = 10**-3
iterations = 10000
result = steepest_descent_constant_step(f, grad, x0, iterations, error_point, error_grad, h)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: constant \ step \ size \ gradient}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
result = steepest_descent_normalized_step(f, grad, x0, iterations, error_point, error_grad, h)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: normalized \ step \ size \ gradient}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Règle d'Armijo et règle de Wolfe¶
Maintenant, on définit la règle d'Armijo et la règle de Wolfe. On introduit également les descentes de gradient (steepest descent gradient en anglais) associées. On rappelle que les règles d'Armijo et de Wolfe permettent de donner une approximation très grossière de la recherche linéaire exacte qui donne lieu à la descente de gradient à pas optimal.
On rappelle la règle d'Armijo.
Soit $\beta \in ]0,1[$, $c \in ]0, \frac{1}{2}[$ et $L \in \mathbb{R}_+$ la constante de Lipschitz de $J$ restreinte à $K$ un compact.
Proposition $\rho_t = -\frac{\langle \nabla J(x^k),d^k \rangle}{L \|d^k\|^2}$ ,
Test d'acceptation Si $J(x^k + \rho_t d^k) \leq J(x^k) + c\, \rho_t\, \langle \nabla J(x^k),d^k \rangle, $ alors $\rho^k = \rho_t$ sinon on reprend le test avec $\beta \rho_t$.
On rappelle la règle de Wolfe.
Soit $c_1 \in ]0, \frac{1}{2}[$, $c_2 \in ]0,\frac{1}{2}[$, et $L \in \mathbb{R}_+$ la constante de Lipschitz de $J$ restreinte à $K$ un compact. Soit $(m_0,M_0) \in \mathbb{R}_+^2$. On fixe également un nombre maximal d'itérations dans la règle de Wolfe ($10^3$ dans notre cas).
Proposition $\alpha_t = -\frac{\langle \nabla J(x^k),d^k \rangle}{L \|d^k\|^2}$,
Test d'acceptation (1) Si $J(x^k + \alpha_t d^k) \leq J(x^k) + m\, \alpha_t\, \langle \nabla J(x^k),d^k \rangle,$ et également ${d^k}^T \nabla J(x^k + \alpha_t d^k) \le -c_2 {d^k}^T \nabla J(x^k)$ alors $\alpha^k = \alpha_t$,
Test d'acceptation (2) Si $J(x^k + \alpha_t d^k) > J(x^k) + m\, \alpha_t\, \langle \nabla J(x^k),d^k \rangle,$ alors on pose $\alpha_t = \frac{m_0+M_0}{2}$ ainsi que $M_0 = \alpha_t$ et on reprend le premier test d'acceptation,
Test d'acceptation (3) Si aucune de ces conditions n'est remplie alors on peut poser $\alpha_t = \frac{m_0+M_0}{2}$ et $m_0 = \alpha_t$ et on reprend le premier test d'acceptation.
Implémentation des règles et descentes de gradient associées¶
In [12]:
def armijo_rule(alpha_0,x,f,f_x,grad_x,d_x,c,beta): #d_x est la direction de descente d_x . grad_x <= 0
# test f(x_new) \leq f(x_0) + c alpha ps{d_x}{grad_x}
test = 1
alpha = alpha_0
while test:
x_new = x+alpha*d_x;
if (f(x_new)<=f_x+c*alpha*np.dot(grad_x,d_x)):
test = 0
else:
alpha = alpha*beta
return alpha
def steepest_descent_armijo(f, grad, x0, iterations, error_point, error_grad, c=0.1,L=100,beta=0.5):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = x0
x_old = x
grad_x = grad(x)
d_x = -grad_x
f_x = f(x)
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = armijo_rule(alpha_0,x,f,f_x,grad_x,d_x,c,beta)
for i in xrange(iterations):
x = x + h*d_x
grad_x = grad(x)
f_x = f(x)
d_x = -grad_x
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = armijo_rule(alpha_0,x,f,f_x,grad_x,d_x,c,beta)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
def wolfe_rule(alpha_0,x,f,grad,f_x,grad_x,d_x,c_1,c_2): #d_x est la direction de descente d_x . grad_x <= 0
# test f(x_new) \leq f(x_0) + c_1 alpha ps{d_x}{grad_x} et \ps{x_new}{d_x} \geq c_2 \ps{x_0}{d_x}
# sinon alpha <- alpha * beta
# On cherche au fur et mesure un opt dans [minorant, majorant]
test = 1
iteration = 0
alpha = alpha_0
minorant = 0
majorant = 1000
while (test)&(iteration<=1000):
x_new = x+alpha*d_x;
if (f(x_new)<=f_x+c_1*alpha*np.dot(grad_x,d_x))&(np.dot(grad(x_new),d_x) >= c_2*np.dot(grad_x,d_x) ):
test = 0
elif (f(x_new)>f_x+c_1*alpha*np.dot(grad_x,d_x)):
majorant = alpha
alpha = (majorant + minorant)/2
iteration = iteration +1
else:
minorant = alpha
alpha = (majorant + minorant)/2
iteration = iteration +1
return alpha
def steepest_descent_wolfe(f, grad, x0, iterations, error_point, error_grad, c_1=0.1,c_2=0.9,L=100):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = x0
x_old = x
grad_x = grad(x)
d_x = -grad_x
f_x = f(x)
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = wolfe_rule(alpha_0,x,f,grad,f_x,grad_x,d_x,c_1,c_2)
for i in xrange(iterations):
x = x + h*d_x
grad_x = grad(x)
f_x = f(x)
d_x = -grad_x
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = wolfe_rule(alpha_0,x,f,grad,f_x,grad_x,d_x,c_1,c_2)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
Application à la fonction de Rosenbrock¶
On teste maintenant les deux règles sur le problème de Rosenbrock.
In [13]:
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
h = 10**-2
iterations = 10000
result = steepest_descent_armijo(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: steepest \ + \ Armijo \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
h = 10**-2
iterations = 10000
result = steepest_descent_wolfe(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: steepest \ + Wolfe \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Application à une fonction quadratique¶
De la même manière, on teste les deux règles et les algorithmes associés sur la fonction quadratique introduite dans la première partie.
In [14]:
M = 5.
m = 1.
n = 3
(f,grad,hess) = mk_quad(M,m,n)
x0 = np.array([5.,5.,5.])
error_point = 10**-10
error_grad = 10**-10
h = 10**-3
iterations = 10000
result = steepest_descent_armijo(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: steepest \ descent \ + \ Armijo \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
result = steepest_descent_wolfe(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: steepest \ descent \ + \ Wolfe \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Gradient conjugué¶
On rappelle ici l'heuristique du gradient conjugué. Si la convergence en temps fini est assurée pour les problèmes quadratiques, ce résultat n'est plus valable sur une fonctionnelle quelconque. On présente ici la version de Polak-Ribière
- Direction de descente : $\beta^{k} = \frac{\langle \nabla J(x^k) - \nabla J(x^{k-1}), \nabla J(x^k) \rangle}{\| \nabla J(x^{k-1}) \|^2}$,
- Pas de descente : on applique une règle de \textit{recherche linéaire inexacte} (ici la règle d'Armijo ou la règle de Wolfe) pour trouver le meilleur pas de descente pour la direction de descente $d^k = -\nabla J(x^k) + \beta^k d^{k-1}$,
- Itération : $x^{k+1} = x^k + \alpha^k d^k$.
Implémentation des méthodes¶
In [15]:
def conjugate_gradient_armijo(f, grad, x0, iterations, error_point, error_grad, c=0.1,L=100,beta=0.5):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = x0
x_old = x
grad_x = grad(x)
d_x = -grad_x
f_x = f(x)
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = armijo_rule(alpha_0,x,f,f_x,grad_x,d_x,c,beta)
for i in xrange(iterations):
x = x + h*d_x
grad_x_old = grad_x
grad_x = grad(x)
f_x = f(x)
kappa = np.dot(grad_x - grad_x_old, grad_x)/np.power(np.linalg.norm(grad_x),2)
d_x = kappa*d_x -grad_x
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = armijo_rule(alpha_0,x,f,f_x,grad_x,d_x,c,beta)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
def conjugate_gradient_wolfe(f, grad, x0, iterations, error_point, error_grad, c_1=0.1,c_2=0.9,L=100):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = x0
x_old = x
grad_x = grad(x)
d_x = -grad_x
f_x = f(x)
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = wolfe_rule(alpha_0,x,f,grad,f_x,grad_x,d_x,c_1,c_2)
for i in xrange(iterations):
x = x + h*d_x
grad_x_old = grad_x
grad_x = grad(x)
f_x = f(x)
kappa = np.dot(grad_x - grad_x_old, grad_x)/np.power(np.linalg.norm(grad_x),2)
d_x = kappa*d_x -grad_x
alpha_0 = -(1./L)*np.dot(d_x,grad_x)/np.power(np.linalg.norm(d_x),2)
h = wolfe_rule(alpha_0,x,f,grad,f_x,grad_x,d_x,c_1,c_2)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}, iteration={}, f(x)={}".format(error_point_list[i], error_grad_list[i],i+1,f(x))
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
Application à la fonction de Rosenbrock¶
On teste maintenant les méthodes de gradient conjugué avec les deux règles introduites sur la fonction de Rosenbrock.
In [16]:
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
iterations = 10000
result = conjugate_gradient_armijo(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: conjugate \ gradient \ + \ Armijo \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
f = rosenbrock
grad = rosenbrock_grad
x0 = np.array([1.1,2.1])
error_point = 10**-10
error_grad = 10**-10
iterations = 10000
result = conjugate_gradient_wolfe(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([1],[1],'g+')
plt.title(r'$\mathrm{Rosenbrock \ minimization: conjugate \ gradient \ + \ Wolfe \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Application à une fonction quadratique¶
Encore une fois, on teste nos algorithmes sur la fonction quadratique introduite précédemment.
In [17]:
M = 5.
m = 1.
n = 3
(f,grad,hess) = mk_quad(M,m,n)
x0 = np.array([5.,5.,5.])
error_point = 10**-10
error_grad = 10**-10
iterations = 10000
result = conjugate_gradient_armijo(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: conjugate \ gradient \ and \ Armijo \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
result = conjugate_gradient_wolfe(f, grad, x0, iterations, error_point, error_grad)
x_list = result['x_list']
all_x_i = np.append(x0[0], x_list[0,:])
all_y_i = np.append(x0[1], x_list[1,:])
plt.plot(all_x_i, all_y_i, 'k+-')
plt.plot(x0[0],x0[1],'r+')
plt.plot([0],[0],'g+')
plt.title(r'$\mathrm{Quadratic \ minimization: conjugate \ gradient \ and \ Wolfe \ rule}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
ATTENTION : ici on n'a pas mis à jour les paramètres des différentes règles. Il convient tout de même d'avoir une idée de l'ordre de grandeur de la constante de Lipschitz $L$ si on veut espérer pouvoir converger (ce qui n'est pas le cas dans le cadre du gradient conjugué avec règle de Wolfe). En modifiant $M$ (par exemple $M = 100$) on obtient une bonne estimation de la constante de Lipschitz avec les paramètres par défaut de l'algorithme. Dans ce cas, les deux algorithmes convergent.
Problème de Lennard-Jones¶
On s'intéresse ici au problème de Lennard-Jones qui donne la configuration d'énergie minimal de certains molécules en utilisant le potentiel de van der Walls pour modéliser les intéractions éléctrostatiques. Ici l'idée est simplement de visualiser les différentes méthodes de descente de gradient. Les fonctions étudiées étant non-convexes on n'a pas d'assurance de convergence. Néanmoins, on peut comparer les performances des différents algorithmes.
Le potentiel de van der Walls¶
<a id = 'vdw'></a>
On commence par donner les fonctions dont on aura besoin pour implémenter nos méthodes de gradient.
In [5]:
def V(x):
return x ** -12 - 2 * x ** - 6
def V2(x):
return x ** -6 - 2 * x ** - 3
def V2der(x):
return -6 * x ** -7 + 6 * x ** -4
def J(u):
N = len(u) / 3
u_v = np.reshape(u, (N, 3))
M = np.zeros([N, N, 3])
M -= u_v
M = M - np.transpose(M, (1, 0, 2))
M = np.sum(M ** 2, 2)
np.fill_diagonal(M, 1)
M = V2(M)
np.fill_diagonal(M, 0)
return .5 * np.sum(M)
def grad_J(u):
N = len(u) / 3
u_v = np.reshape(u, (N, 3))
M = np.zeros([N, N, 3])
M -= u_v
M = M - np.transpose(M, (1, 0, 2))
Mnorm = np.sum(M ** 2, 2)
np.fill_diagonal(Mnorm, 1)
Mnorm = V2der(Mnorm)
np.fill_diagonal(Mnorm, 0)
grad = np.reshape(Mnorm, (N**2, 1)) * np.reshape(M, (N ** 2, 3))
grad = np.reshape(grad, (N, N, 3))
grad = np.sum(grad, 1)
return 2 * np.ravel(grad)
On va maintenant visualiser le potentiel de van der Walls pour mettre en évidence ses variations
In [19]:
r = np.linspace(0.85,1.5,1000)
Vr = V(r)
plt.plot(r,Vr)
plt.xlabel('rayon')
plt.ylabel('Force')
plt.title(r'$\mathrm{Potentiel \ de \ van \ der \ Walls}$')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.show()
Application pour $N=4$ et $N=13$¶
<a id = "appli"></a>
In [4]:
N = 4
error_point = 10 ** -10
error_grad = 10 ** -10
x0 = rnd.random(3*N)
x0[0:2] = 0
iterations = 10000
result = steepest_descent_armijo(J, grad_J, x0, iterations, error_point, error_grad)
init_pos = np.reshape(x0, (N,3))
final_pos = np.reshape(result['x_list'][:,-1], (N,3))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(init_pos[:,0], init_pos[:,1], init_pos[:,2], '^')
ax.scatter(final_pos[:,0], final_pos[:,1], final_pos[:,2], '^')
for i in range(4):
ax.plot([init_pos[i,0], final_pos[i,0]],
[init_pos[i,1], final_pos[i,1]],
[init_pos[i,2], final_pos[i,2]], 'g')
for j in range(4):
ax.plot([final_pos[i,0], final_pos[j,0]],
[final_pos[i,1], final_pos[j,1]],
[final_pos[i,2], final_pos[j,2]], 'r')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
print(J(np.ravel(init_pos)))
print(J(np.ravel(final_pos)))
On tente aussi pour le cas où $N = 13$. Dans ce cas, la configuration minimale est un icosaèdre avec son centre de gravité.
In [21]:
N = 13
error_point = 10 ** -10
error_grad = 10 ** -10
x0 = rnd.random(3*N)
x0[0:2] = 0
iterations = 100000
result = steepest_descent_armijo(J, grad_J, x0, iterations, error_point, error_grad)
init_pos = np.reshape(x0, (N,3))
final_pos = np.reshape(result['x_list'][:,-1], (N,3))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(final_pos[:,0], final_pos[:,1], final_pos[:,2], '^')
for i in range(13):
ax.plot([init_pos[i,0], final_pos[i,0]],[init_pos[i,1], final_pos[i,1]],[init_pos[i,2], final_pos[i,2]], 'g')
moy = np.mean(final_pos,0)
d = np.linalg.norm(final_pos - moy,axis = 1)
center = np.argmin(d)
vertices = final_pos[np.arange(len(final_pos))!=center,:]
vertdiff = np.zeros([N-1, N-1, 3])
vertdiff -= vertices
vertdiff = vertdiff - np.transpose(vertdiff, (1, 0, 2))
vertnorm = np.sum(vertdiff ** 2, 2)
for i in range(12):
indsort = np.argsort(vertnorm[:,i])
for j in range(5):
js = indsort[j]
ax.plot([vertices[i,0], vertices[js,0]],[vertices[i,1], vertices[js,1]],[vertices[i,2], vertices[js,2]], 'r')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
print J(np.ravel(init_pos))
print J(np.ravel(final_pos))
ATTENTION : comme tout à l'heure le problème est non-convexe et donc on peut être coincé dans des minimums locaux. Il peut être nécessaire de relancer plusieurs fois l'algorithme avec des initialisations différentes pour obtenir un icosaèdre. La valeur minimale est environ $-44.32$.
Méthode de Levenberg-Marquardt, Gauss-Newton et application à la régression non-linéaire¶
Cette partie du TP est fortement inspirée de l'application proposée dans le Scipy Cookbook.
Préambule : régression non-linéaire¶
Dans cette première partie, on définit les fonctions qui vont nous être utiles pour tester une méthode du second order, la méthode de Gauss-Newton et une méthode de quasi-Newton, la méthode de Levenberg-Marquardt. On commence par introduire les fonctions et générer des données.
In [22]:
def mk_wave(A, sigma, omega):
def wave_fun(x):
return A * np.exp(-sigma * x) * np.sin(omega * x)
def wave_grad(x):
dim = np.max(np.shape(x))
grad = np.zeros([3, dim])
grad[0, :] = np.sin(omega * x)
grad[1, :] = -A * x * np.sin(omega * x)
grad[2, :] = A * x * np.cos(omega * x)
grad = np.exp(-sigma * x) * grad
return grad
def wave_hessian(x):
dim = np.max(np.shape(x))
s = np.sin(omega * x) * x * np.exp(-sigma * x)
c = np.cos(omega * x) * x * np.exp(-sigma * x)
scal = A * x ** 2 * np.exp(-sigma * x)
hessian = np.zeros([3, 3, dim])
hessian[0, 1, :] = -x* s
hessian[1, 0, :] = -x * s
hessian[0, 2, :] = hessian[2, 0, :] = x * c
hessian[1, 1, :] = scal * s
hessian[2, 2, :] = -scal * s
hessian[1, 2, :] = hessian[2, 1, :] = -scal * c
return hessian
return wave_fun, wave_grad, wave_hessian
def generate_data(x, A, sigma, omega, noise=0, n_outliers=0, random_state=0):
y = A * np.exp(-sigma * x) * np.sin(omega * x)
rnd = np.random.RandomState(random_state)
error = noise * rnd.randn(x.size)
outliers = rnd.randint(0, x.size, n_outliers)
error[outliers] *= 35
return y + error
A = 5
sigma = 0.1
omega = 0.1 * 2 * np.pi
param_true = np.array([A, sigma, omega])
wave_fun, wave_grad, wave_hessian = mk_wave(A,sigma,omega)
noise = 0.1
x_min = 0
x_max = 20
x_train = np.linspace(x_min, x_max, 30)
y_train = generate_data(x_train, A, sigma, omega, noise=noise, n_outliers=5)
plt.plot(x_train, y_train, 'o', label = 'data')
plt.plot(x_train, wave_fun(x_train), 'r-', label = 'signal')
plt.title(r'$\mathrm{Data \ with \ outliers}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.legend()
plt.show()
In [23]:
r = np.linspace(0, 5, 100)
linear = r**2
huber = r**2
huber[huber > 1] = 2 * r[huber > 1] - 1
soft_l1 = 2 * (np.sqrt(1 + r**2) - 1)
cauchy = np.log1p(r**2)
arctan = np.arctan(r**2)
plt.plot(r, linear, label='linear')
plt.plot(r, huber, label='Huber')
plt.plot(r, soft_l1, label='soft_l1')
plt.plot(r, cauchy, label='Cauchy')
plt.plot(r, arctan, label='arctan')
plt.xlabel("$r$")
plt.ylabel(r"$\rho(r^2)$")
plt.title(r'$\mathrm{Different \ loss \ functions}$')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.legend()
plt.show()
On en profite pour définir les gradients et hessiennes qui nous seront utiles pour mettre en place la méthode de Newton et la méthode de Levenberg-Marquardt. On se concentre sur la fonction de coût quadratique.
In [6]:
def mk_nonlinreg(x_train, y_train):
def nonlinreg_fun(param):
wave_fun = mk_wave(param[0], param[1], param[2])[0]
return np.sum((wave_fun(x_train) - y_train) ** 2)
def nonlinreg_grad(param):
wave_fun, wave_grad = mk_wave(param[0], param[1], param[2])[:2]
grad = 2 * wave_grad(x_train) * (wave_fun(x_train) - y_train)
return np.sum(grad, 1)
def nonlinreg_hessian(param, method='newton', mu=0.1):
wave_fun, wave_grad, wave_hessian = mk_wave(
param[0], param[1], param[2])
grad = wave_grad(x_train)
if method == 'newton':
hess1 = 2 * wave_hessian(x_train)*\
(wave_fun(x_train) - y_train)
hess1 = np.sum(hess1, 2)
elif method == 'lm':
hess1 = mu * np.identity(3)
hess2 = 2 * np.dot(grad, grad.T)
return hess1 + hess2
return nonlinreg_fun, nonlinreg_grad, nonlinreg_hessian
Implémentation de la méthode de Gauss-Newton et de Levenberg-Marquardt¶
On introduit maintenant la méthode de Newton et la méthode de Levenberg-Marquardt. Il s'agit simplement de modifier la direction de descente pour prendre en compte une information sur l'évolution de la direction de descente (codée par la Hessienne). Un bon moyen de comprendre pourquoi l'introduction de la hessienne est d'introduire les régularisations quadratiques suivantes :
- $f_1(x) = f(x^k) + \nabla f(x^k)(x-x^k) + \frac{1}{2 \alpha}\|x-x^k\|^2$
- $f_2(x) = f(x^k) + \nabla f(x^k)(x-x^k) + \frac{1}{2}(x-x^k)^T H(x^k) (x-x^k)$
La minimisation de la première fonctionnelle donne la méthode de descente de gradient à pas fixe tandis que la seconde donne la méthode de Newton. Néanmoins on a besoin du fait que $H(x^k)$ soit une forme quadratique positive, ce qui est acquis si la fonction $f$ est convexe.
In [25]:
def newton_descent(f, grad, hessian, x0, iterations, error_point, error_grad, h, method = 'newton', mu = 1):
dim = np.max(np.shape(x0))
x_list = np.zeros([dim,iterations])
f_list = np.zeros(iterations)
error_point_list = np.zeros(iterations)
error_grad_list = np.zeros(iterations)
x = np.asarray(x0)
x_old = x
grad_x = grad(x)
count = 0
for i in xrange(iterations):
M = hessian(x,method,mu)
if any(v<0 for v in np.linalg.eig(M)[0]):
count += 1
descent = np.linalg.solve(M, grad(x))
x = x - h*descent
grad_x = grad(x)
f_x = f(x)
x_list[:,i] = x
f_list[i] = f_x
error_point_list[i] = np.linalg.norm(x - x_old)
error_grad_list[i] = np.linalg.norm(grad_x)
if i % 1000 == 0:
print "iter={}, x={}, f(x)={}".format(i+1, x, f(x))
if (error_point_list[i] < error_point)|(error_grad_list[i] < error_grad):
break
x_old = x
print "point error={}, grad error={}".format(error_point_list[i], error_grad_list[i])
print "number of negative eigenvalues={}".format(count)
return { 'x_list' : x_list[:,0:i], 'f_list' : f_list[0:i], 'error_point_list' : error_point_list[0:i], 'error_point_list' : error_point_list[0:i]}
Dans notre application on compare les algorithmes de Newton ou quasi-Newton aux algorithmes de descente de gradient introduits plutôt.
In [26]:
f, grad, hess = mk_nonlinreg(x_train, y_train)
x0 = rnd.random(3)
true_hessian = hess(x0)
lm_hessian = hess(x0, 'lm', 0.1)
print "true cond = {}, lm cond = {}".format(np.linalg.cond(true_hessian), np.linalg.cond(lm_hessian))
error_point = 10 ** -10
error_grad = 10 ** -10
h = 10 ** -3
iterations = 10000
result_newton = newton_descent(f, grad, hess, x0, iterations, error_point, error_grad, h)
result_lm = newton_descent(f, grad, hess, x0, iterations, error_point, error_grad, h, method = 'lm', mu = 10 ** -1)
result_wolfe = steepest_descent_wolfe(f, grad, x0, iterations, error_point, error_grad)
(param_newton, param_lm, param_wolfe) = (result['x_list'][:,-1] for result in (result_newton, result_lm, result_wolfe))
print "Newton parameters = {}, Levenberg-Marquardt parameters = {}, Wolfe parameters = {}".format(param_newton, param_lm, param_wolfe)
newton_wave_fun = mk_wave(param_newton[0], param_newton[1], param_newton[2])[0]
lm_wave_fun = mk_wave(param_lm[0], param_lm[1], param_lm[2])[0]
wolfe_wave_fun = mk_wave(param_wolfe[0], param_wolfe[1], param_wolfe[2])[0]
plt.plot(x_train, y_train, 'o', label = 'data')
plt.plot(x_train, wave_fun(x_train), 'r-', label = 'signal')
plt.plot(x_train, newton_wave_fun(x_train), 'c-', label = 'Newton')
plt.plot(x_train, lm_wave_fun(x_train), 'g-', label = 'Levenberg')
plt.plot(x_train, wolfe_wave_fun(x_train), 'b-', label = 'Wolfe')
plt.title(r'$\mathrm{Data \ with \ outliers}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.legend()
plt.show()
On constate que la méthode de Newton échoue. En effet, la forme quadratique mise en jeu n'est pas positive.
Complément : test sur le problème de régression linéaire¶
En complément, on montre que la méthode de Gauss-Newton converge bien dans le cas où la fonction que l'on cherche à interpoler est convexe en les paramètres. C'est le cas de la fonction linéaire $f: \ x \mapsto ax+b$ où les paramètres sont $a$ et $b$. On reprend la même étude que précédemment.
In [27]:
def mk_lin(a, b):
def lin_fun(x):
return a * x + b
def lin_grad(x):
grad = np.ones([2, x.size])
grad[0,:] = x
return grad
def lin_hessian(x, method = 'newton', mu = 0.1):
return np.zeros([2,2,x.size])
return lin_fun, lin_grad, lin_hessian
def mk_linreg(x_train, y_train):
def linreg_fun(param):
lin_fun = mk_lin(param[0], param[1])[0]
return np.sum((lin_fun(x_train) - y_train) ** 2)
def linreg_grad(param):
lin_fun, lin_grad = mk_lin(param[0], param[1])[:2]
grad = 2 * lin_grad(x_train) * (lin_fun(x_train) - y_train)
return np.sum(grad, 1)
def linreg_hessian(param, method='newton', mu=0.1):
lin_fun, lin_grad, lin_hessian = mk_lin(
param[0], param[1])
grad = lin_grad(x_train)
if method == 'newton':
hess1 = 2 * lin_hessian(x_train) * (lin_fun(x_train) - y_train)
hess1 = np.sum(hess1, 2)
elif method == 'lm':
hess1 = mu * np.identity(2)
hess2 = 2 * np.dot(grad, grad.T)
return hess1 + hess2
return linreg_fun, linreg_grad, linreg_hessian
In [28]:
a = 2
b = 3
lin_fun, lin_grad, lin_hessian = mk_lin(a,b)
noise = 5
x_min = 0
x_max = 20
n_data = 50
x_train = np.linspace(x_min, x_max, n_data)
y_train = a * x_train + b + noise * rnd.randn(x_train.size)
plt.plot(x_train, y_train, 'o', label = 'data')
plt.plot(x_train, lin_fun(x_train), 'r-', label = 'signal')
plt.title(r'$\mathrm{Data \ with \ outliers}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.legend()
plt.show()
f, grad, hess = mk_linreg(x_train, y_train)
error_point = 10 ** -10
error_grad = 10 ** -10
h = 10 ** -2
iterations = 1000
x0 = rnd.random(2)
result_newton = newton_descent(f, grad, hess, x0, iterations, error_point, error_grad, h)
result_lm = newton_descent(f, grad, hess, x0, iterations, error_point, error_grad, h, method = 'lm', mu = 0.01)
result_wolfe = steepest_descent_wolfe(f, grad, x0, iterations, error_point, error_grad)
(param_newton, param_lm, param_wolfe) = (result['x_list'][:,-1] for result in (result_newton, result_lm, result_wolfe))
print "Newton parameters = {}, Levenberg-Marquardt parameters = {}, Wolfe parameters = {}".format(param_newton, param_lm, param_wolfe)
newton_lin_fun, newton_lin_grad, newton_lin_hessian = mk_lin(param_newton[0], param_newton[1])
lm_lin_fun, lm_lin_grad, lm_lin_hessian = mk_lin(param_lm[0], param_lm[1])
wolfe_lin_fun, wolfe_lin_grad, wolfe_lin_hessian = mk_lin(param_wolfe[0], param_wolfe[1])
plt.plot(x_train, y_train, 'o', label = 'data')
plt.plot(x_train, lin_fun(x_train), 'r-', label = 'signal')
plt.plot(x_train, newton_lin_fun(x_train), 'c-', label = 'Newton')
plt.plot(x_train, lm_lin_fun(x_train), 'g-', label = 'Levenberg')
plt.plot(x_train, wolfe_lin_fun(x_train), 'b-', label = 'Wolfe')
plt.title(r'$\mathrm{Data \ with \ outliers}$')
plt.xlabel('x')
plt.ylabel('y')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.legend()
plt.show()
