opti-non-lin/cours3.ipynb

In [138]:

import numpy as np
import matplotlib.pyplot as plt
import time

Méthode Newton $\mathbb{R}$

In [139]:

def f1(x):
    f = 2*x-np.log(x)
    g = 2 - 1/x
    h = 1/(x**2)
    return f, g, h

In [140]:

def f2(x):
    f = x**3 - 6*x**2 + 4*x + 2
    g = 3*x**2 - 12*x + 4
    h = 6*x - 12
    return f, g, h

In [141]:

def NewtonR(fct, x0, maxiter, tol=1e-9):
    x = np.zeros(maxiter)
    x[0] = x0

    for k in range(maxiter-1):
        f, g, h = fct(x[k])
        x[k+1] = x[k] - g/h

        if abs(g) < tol:
            return x[:k+1]

    return x

In [142]:

# Ici ça marche super bien

absc = np.linspace(0.01, 1.5, 1000)
oord = f1(absc)[0]

x0 = 0.1

x = NewtonR(f1, x0, 1000)
for k in range(len(x)):
    f, g, h = f1(x[k])
    plt.plot(x[k], f, 'ob')

plt.plot(absc, oord, 'g')

Out[142]:

[<matplotlib.lines.Line2D at 0x1eec176bb10>]

In [143]:

# Ici ça dépend du point initial

absc = np.linspace(-2, 5, 1000)
oord = f2(absc)[0]

x0 = [3,0]
c = ['b', 'r']

for i, x0 in enumerate(x0):
    x = NewtonR(f2, x0, 1000)
    for k in range(len(x)):
        f, g, h = f2(x[k])
        plt.plot(x[k], f, 'o' + c[i])

plt.plot(absc, oord, 'g')

Out[143]:

[<matplotlib.lines.Line2D at 0x1eec408c050>]

Examen d'implémentation | Gradient et Newton

$f(x_1,x_2)=2x_1^2-1.05x_1^4+x_1^6/6+x_1x_2+x_2^2$

$\nabla f(x_1,x_2)= \begin{pmatrix} 4x_1 - 4.2x_1^3 + x_1^5 + x_2 & x_1 + 2x_2 \\ \end{pmatrix}$

$\nabla^2 f(x_1,x_2)= \begin{pmatrix} 4 - 12.6x_1^2 + 5x_1^4 & 1 \\ 1 & 2 \\ \end{pmatrix}$

In [144]:

def fexamen(x, nargout=2):
    x1 = x[0]
    x2 = x[1]
    f = 2*x1**2 - 1.05*x1**4 + x1**6/6 + x1*x2 + x2**2 

    g = np.array([4*x1 - 4*1.05*x1**3 + x1**5 + x2, x1 + 2*x2])

    if nargout==3:
        H = np.array([[4 - 12*1.05*x1**2 + 5*x1**4, 1],
                      [              1            , 2]])
        return f, g, H

    return f, g

In [145]:

# Vérification de fexamen: OK
x = np.array([0,1])
print(fexamen(x))

(np.float64(1.0), array([1., 2.]))

In [146]:

def w1(f,x,d,alpha,beta1):
    f0, g0 = f(x) # à l'itéré courant
    f1, g1 = f(x+alpha*d) # à l'itéré potentiel
    return f1 <= f0 + alpha*beta1 * np.dot(d, g0)

def w2(f,x,d,alpha,beta2):
    f0, g0 = f(x) # à l'itéré courant
    f1, g1 = f(x+alpha*d) # à l'itéré potentiel   
    return np.dot(d, g1) >= beta2*np.dot(d, g0)

def wolfebissection(f,x,d,alpha,beta1=1e-3,beta2=0.9):
    aleft = 0
    aright = np.inf

    while True:
        w_1 = w1(f,x,d,alpha,beta1)
        w_2 = w2(f,x,d,alpha,beta2)

        if w_1 and w_2:
            break

        elif not w_1:
            aright = alpha
            alpha = (aleft+aright)/2

        elif not w_2:
            aleft = alpha
            if aright<np.inf:
                alpha = (aleft+aright)/2
            else:
                alpha *= 2 # 2 est arbitraire

    return alpha

In [147]:

def methgradient(f,x0,maxiter=300,tol=1e-6):
    x = x0.copy()
    fobj = np.zeros(maxiter)
    ngrad = np.zeros(maxiter)

    for k in range(maxiter):
        fx, g = f(x)
        d = -g
        alpha = wolfebissection(f,x,d,1)

        x = x+alpha*d # calcul nouvel itéré

        # sauvegarde des f et g
        fx, g = f(x)
        fobj[k] = fx
        ngrad[k] = np.linalg.norm(g)

        # critère d'arrêt
        if np.linalg.norm(g) <= tol:
            break

    return x, fobj[:k+1], ngrad[:k+1]

In [148]:

def methNewton(f,x0,maxiter=300,tol=1e-6):
    x = x0.copy()
    fobj = np.zeros(maxiter)
    ngrad = np.zeros(maxiter)

    for k in range(maxiter):
        fx, g, H = f(x, nargout=3) # besoin de H
        z = np.linalg.solve(H,g)

        x = x-z # calcul nouvel itéré

        # sauvegarde des f et g
        fx, g = f(x)
        fobj[k] = fx
        ngrad[k] = np.linalg.norm(g)

        # critère d'arrêt
        if np.linalg.norm(g) <= tol:
            break

    return x, fobj[:k+1], ngrad[:k+1]

In [149]:

x0 = np.array([0.1, 0.1])

xg, fg, ngg = methgradient(fexamen,x0)
xn, fn, ngn = methNewton(fexamen,x0)

plt.semilogy(fg, label="gradient")
plt.semilogy(fn, label="Newton")
plt.title("Evolution de la fonction objectif")
plt.legend()
plt.show()

plt.semilogy(ngg, label="gradient")
plt.semilogy(ngn, label="Newton")
plt.title("Evolution de la norme du gradient")
plt.legend()
plt.show()

In [150]:

print(f"Avec la méthode du gradient : x={xg}")
f, g, H = fexamen(xg, nargout=3)
print(f"fobj={f} | grad={g} | H={H}")
print(np.linalg.eig(H))
# On a un point ou le gradient vaut 0 et la matrice hessienne est definie positive
# => on a un minimum local

print(f"Avec la méthode de Newton : x={xn}")
f, g, H = fexamen(xn, nargout=3)
print(f"fobj={f} | grad={g} | H={H}")
print(np.linalg.eig(H))
# On a un point ou le gradient vaut 0 et la matrice hessienne est indéfinie
# => on a un point selle

Avec la méthode du gradient : x=[-9.13230123e-08  9.55910628e-08]
fobj=1.7087772621488852e-14 | grad=[-2.69700986e-07  9.98591133e-08] | H=[[4. 1.]
 [1. 2.]]
EigResult(eigenvalues=array([4.41421356, 1.58578644]), eigenvectors=array([[ 0.92387953, -0.38268343],
       [ 0.38268343,  0.92387953]]))
Avec la méthode de Newton : x=[ 3.64930450e-08 -1.82465225e-08]
fobj=2.3305490890415854e-15 | grad=[1.27725658e-07 4.33680869e-19] | H=[[4. 1.]
 [1. 2.]]
EigResult(eigenvalues=array([4.41421356, 1.58578644]), eigenvectors=array([[ 0.92387953, -0.38268343],
       [ 0.38268343,  0.92387953]]))

$f(x)=(x_1-2x_2)^4 + (x_2-3x_3)^4 + (x_3-4x_4)^4 + \ldots + (x_n-(x_1+1))^4$

$x_1 = 2x_2 = 2*3*4(x_1+1) => -23x_1=24 => x_1=-24/23$

$x_2 = 3x_3 = 3*4(x_1+1) => x_2 = \ldots$

$\ldots$

$x_i = \frac{-n!}{i(n!-1)}$

In [151]:

def fsize(x,nargout=2):
  n = len(x)

  #Calcul de f
  v      = 1.0*x
  v[:-1] = v[:-1]-x[1:]*np.arange(2,1+len(x))
  v[-1]  = x[-1]-x[0]-1
  f      = np.sum(v**4)

  #Calcul du gradient uniquement s'il est demandé
  g     = 4*v**3
  g[1:] = g[1:]-4*(v[:-1]**3)*np.arange(2,1+len(x))
  g[0]  = g[0]-4*v[-1]**3

  #Calcul de la hessienne uniquement si elle est demandée
  if nargout==3:
    H = np.zeros((n,n))
    for i in range(1,n):
      H[i-1:i+1,i-1:i+1] = H[i-1:i+1,i-1:i+1]+ 12*v[i-1]**2*np.array([[1,-(i+1)],[-(i+1),(i+1)**2]])
    H[0,0]   += 12*v[-1]**2
    H[0,-1]  -= 12*v[-1]**2
    H[-1,0]  -= 12*v[-1]**2
    H[-1,-1] += 12*v[-1]**2
  if nargout==3:
    return f, g, H
  else:
    return f, g

In [152]:

def methgradientTime(f,x0,maxtime,maxiter=300):
    start_time = time.time()
    x = x0.copy()
    fobj = np.zeros(maxiter)
    temps = np.zeros(maxiter)

    for k in range(maxiter):
        fx, g = f(x)
        d = -g
        alpha = wolfebissection(f,x,d,1)

        x = x+alpha*d # calcul nouvel itéré

        # sauvegarde des f et g
        fx, g = f(x)
        fobj[k] = fx
        temps[k] = time.time() - start_time

        if temps[k] > maxtime:
            break

        if np.linalg.norm(g) <= 1e-6:
            break

    return x, fobj[:k+1], temps[:k+1]

In [153]:

def methNewtonTime(f,x0,maxtime,maxiter=300):
    start_time = time.time()
    x = x0.copy()
    fobj = np.zeros(maxiter)
    temps = np.zeros(maxiter)

    for k in range(maxiter):
        fx, g, H = f(x, nargout=3) # besoin de H
        z = np.linalg.solve(H,g)

        x = x-z # calcul nouvel itéré

        # sauvegarde des f et g
        fx, g = f(x)
        fobj[k] = fx
        temps[k] = time.time() - start_time

        if temps[k] > maxtime:
            break

        # critère d'arrêt
        if np.linalg.norm(g) <= 1e-6:
            break

    return x, fobj[:k+1], temps[:k+1]

In [ ]:

# La méthode de Newton a le temps de calculer moins d'itérés que le gradient
# à cause du coût de calcul plus élevé. En terme de temps, le gradient est meilleur.

n = 5000
x0 = np.ones(n)

maxtime = 15
xg, fg, tg = methgradientTime(fsize, x0, maxtime)
xn, fn, tn = methNewtonTime(fsize, x0, maxtime)

plt.semilogy(fg, label="gradient")
plt.semilogy(fn, label="Newton")
plt.title("Evolution de la fonction objectif en fonction des itérations")
plt.legend()
plt.show()

plt.semilogy(tg, fg, label="gradient")
plt.semilogy(tn, fn, label="Newton")
plt.title("Evolution de la norme du gradient, en fonction du temps")
plt.legend()
plt.show()