In [26]:
import numpy as np
import matplotlib.pyplot as plt

def log_1(t):
    if t <= 0:
        return -float("inf")
    return np.log(t)
log = np.vectorize(log_1)

# comme en TP, on fournit des fonctions permettant de vérifier le gradient et/ou la hessienne
def verifier_gradient(f,gradf,x0):
    N = len(x0)
    gg = np.zeros(N)
    for i in range(N):
        eps = 1e-5
        e = np.zeros(N)
        e[i] = eps
        gg[i] = (f(x0+e) - f(x0-e))/(2*eps)
    print('erreur numérique dans le calcul du gradient: %g (doit être petit)' % np.linalg.norm(gradf(x0)-gg))
    
def verifier_hessienne(gradf,hessf,x0):
    N = len(x0)
    H = np.zeros((N,N))
    for i in range(N):
        eps = 1e-5
        e = np.zeros(N)
        e[i] = eps
        H[i,:] = (gradf(x0+e) - gradf(x0-e))/(2*eps)
    print('erreur numerique dans le calcul de la hessienne: %g (doit etre petit)' % np.sum((H-hessf(x0)))**2)
    
# on donne également la fonction pas_armijo/gradient_armijo
def pas_armijo(f,gradf,x,v):
    t = 1
    m = np.dot(v,gradf(x))
    alpha=0.3
    beta=0.5
    while f(x+t*v) > f(x) + alpha*t*m:
        t = beta*t
    return t

def gradient_armijo(f,gradf,x0,err=1e-4):
    # ne pas hésiter à copier-coller et adapter gradient_optimal...
    x = x0
    G = []
    k = 0 # nombre d'itérations
    maxiter=200
    while(True): 
        k = k+1
        if k > maxiter:
            print('erreur: nombre maximum d\'itérations atteint')
            break
        # calcul de d, t, 
        d = -gradf(x)
        G.append(np.linalg.norm(d))
        if np.linalg.norm(d) <= err:
            break
        t = pas_armijo(f,gradf,x,d)
        x = x + t*d
    return x,np.array(G)

Partie 1.2 -- Étude numérique

$$\newcommand{\eps}{\varepsilon}$$

QN1. Écrire deux fonctions f(x), gradf(x) calculant respectivement $f(x),\nabla f(x)$. (On pourra utiliser la fonction verifier_gradient() pour vérifier le calcul)

**(Attention, il faut utiliser la fonction log donnée ci-dessus plutôt que np.log)**

In [27]:
#<
eps = 1e-3
def f(x):
    return .5*np.linalg.norm(x - z)**2 - eps*(log(x[0]) + log(x[1]) + log(1-x[0]-x[1]))
def gradf(x):
    return x - z - eps*(np.array([1/x[0], 1/x[1]]) + 1/(x[0]+x[1]-1)*np.ones(2))
#>
    
verifier_gradient(f,gradf,np.random.rand(2)/3)
erreur numérique dans le calcul du gradient: 3.98877e-10 (doit être petit)

QN2. La fonction x,G = gradient_armijo(f,gradf,x0) donnée en préambule implémente la descente de gradient avec rebroussement d'Armijo. Elle retourne la solution calculée x et un tableau contenant $(\|{d^{(k)}}\|)$. En utilisant cette fonction, tracer sur un même figure et en échelle logarithmique $k\mapsto \|{d^{(k)}}\|$ en partant de $x^{(0)} = (\frac{1}{4},\frac{1}{4})$, pour plusieurs choix du paramètre $\eps \in \{.1,.01,.001,.0001\}$.

In [28]:
#<
for eps in [.1,.01,0.001,1e-4]:
    x,G = gradient_armijo(f,gradf,.25*np.ones(2))
    plt.semilogy(G,label='eps=%g' % eps)
plt.legend()
#>
Out[28]:
<matplotlib.legend.Legend at 0x7fccb9d57cf8>

QN3 Écrire une fonction hessf(x) calculant $\mathrm{D}^2 f(x)$. (On pourra utiliser verifier_hessienne() pour vérifier le calcul)

In [29]:
#<
def hessf(x):
    return np.eye(2) + eps*(np.diag([1/x[0]**2, 1/x[1]**2]) + 1/((1-x[0]-x[1])**2)*np.ones((2,2)))
#>
    
verifier_hessienne(gradf,hessf,np.random.rand(2)/3)
erreur numerique dans le calcul de la hessienne: 1.05158e-21 (doit etre petit)

QN4. Implémenter la méthode de Newton avec rebroussement d'Armijo (appelée (Newton) dans le sujet). La fonction Python x,G = newton_armijo(f,gradf,x0) retournera la solution calculée x et un tableau contenant $(\|{d^{(k)}}\|)$. En utilisant cette fonction, tracer en échelle logarithmique $k\mapsto \|{d^{(k)}}\|$ en partant de $x^{(0)} = (\frac{1}{4},\frac{1}{4})$ et pour plusieurs choix du paramètre $\eps \in \{.1,.01,.001,.0001\}$.

In [30]:
#<
def newton_armijo(f,gradf,hessf,x0,err=1e-4,maxiter=100):
    # ne pas hésiter à copier-coller et adapter gradient_optimal...
    x = x0
    G = []
    k = 0 # nombre d'itérations
    while(True): 
        k = k+1
        if k > maxiter:
            print('erreur: nombre maximum d\'itérations atteint')
            break
        # calcul de d, t, 
        g = gradf(x)
        d = -np.linalg.solve(hessf(x),g)
        # on pourrait aussi utiliser np.linalg.inv (mais ça serait moins efficace)
        # d = -np.dot(np.linalg.inv(hessf(x)),gradf(x))
        G.append(np.linalg.norm(g))
        if np.linalg.norm(g) <= err:
            break
        t = pas_armijo(f,gradf,x,d)
        x = x + t*d
    return x,np.array(G)

for eps in [.1,.01,.001,1e-4]:
    x,G = newton_armijo(f,gradf,hessf,.25*np.ones(2))
    plt.semilogy(G,label='eps=%g' % eps)
plt.legend()
#>
Out[30]:
<matplotlib.legend.Legend at 0x7fccb94f15f8>
In [ ]: