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

# comme en TP, on fournit des fonctions permettant de vérifier 
# les calculs de gradient et/ou de hessienne
def verifier_derivee(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))
    
# effectue la descente de gradient avec pas d'Armijo:
# Arguments: f = fonction à minimiser
#            gradf = fonction calculant le gradient de d
#            x0 = point de départ
# Retourne: x, F où
#           x est la solution trouvée
#           F = [f(x^(0)), f(x^(1)), ...]
def gradient_armijo(f,gradf,x0):
    niter = 200
    x = x0.copy()
    F = np.zeros(niter)
    for k in range(niter):
        d = -gradf(x)
        m = -np.dot(d,d)
        t = 1
        while f(x+t*d) > f(x) + 0.3*t*m:
            t = 0.5*t
        F[k] = f(x)
        x = x + t*d
    return x,F

Partie I: Méthode de Newton pure en dimension $1$

Dans cette partie, on considère la fonction $g$ donnée par la formule (3), i.e.

$$ g(t) = a t^2 + b t + \sqrt{\varepsilon + (c+t)^2} \hbox{ où } \varepsilon>0, a\geq 0 $$

QN1 Écrire trois fonctions g, gprime et gseconde prenant en argument a,b,c,eps,t et retournant respectivement $g(t)$, $g'(t)$ et $g''(t)$.

In [ ]:
def g(a,b,c,eps,t):
    # à compléter
# définir gprime et gseconde


## vérification du calcul de gprime et gseconde
a=0.1; b=2; c=-3; eps=0.01
verifier_derivee(lambda t: g(a,b,c,eps,t),lambda t: gprime(a,b,c,eps,t),np.random.rand(1))
verifier_derivee(lambda t: gprime(a,b,c,eps,t),lambda t: gseconde(a,b,c,eps,t),np.random.rand(1))

QN2 Écrire une fonction newton_pur(a,b,c,eps) qui implémente la méthode de Newton pure donnée par (2): $$ \begin{cases} t^{(0)} = 0,\\ t^{(k+1)} = t^{(k)} - g'(t^{(k)})/g''(t^{(k)}). \end{cases}$$ et retourne la valeur après cinq itérations, i.e. $t^{(5)}$. On utilisera les fonctions gprime(a,b,c,eps,t) pour calculer $g'(t)$.

In [ ]:
def newton_pur(a, b, c, eps):
    # à compléter

## vérification du fonctionnement de la fonction Newton pur
print("cas 1: solution trouvée %g, solution attendue -5.0003905411" % newton_pur(0.1,2,-3,0.01)) 
print("cas 2: solution trouvée %g, solution attendue -995.0000000024029" % newton_pur(0.1,200,-25,0.001))

Partie II Application en dimension $d\geq 1$: problème LASSO

Il s'agit de résoudre le problème de minimisation (4) $$ \min_{x\in \mathbb{R}^d} f(x) \hbox{ où } f(x) = \frac{\mu}{2} \|Ax - y\|^2 + \sum_{1\leq i\leq d} s(x_i) \hbox{ et } s(t) = \sqrt{\varepsilon + t^2}.$$ où $A,y,\mu$ et $\varepsilon$ sont les suivants. Pour vérifier votre code, le minimiseur du problème est stocké dans xsol et le minimum dans fmin

In [ ]:
# Données du problème
A = np.array([[1,0,1,0,1,0,0,0],
              [1,0,0,1,0,1,0,0],
              [0,1,1,0,0,0,1,0],
              [0,1,0,1,1,0,0,1]])
y = np.array([-1,1,-4,2])
n,d = A.shape
eps = 1e-4
mu = 100

# solution du problème, pour vérification!
fmin = 6.02030542092432
xsol = [0.0000017522, -0.000002031,  -1.5035155249, 0.9999950239,  
        0.5035136323, 0.0000016119, -2.486482525,   0.4864954868]

QN3 Écrire deux fonctions f, gradf prenant en argument x et retournant respectivement $f(x)$ et $\nabla f(x)$.

In [ ]:
def f(x):
    # à compléter
# définir gradf

## vérification du gradient de f
verifier_derivee(f,gradf,np.random.rand(d))

QN4 Utiliser la fonction gradient_armijo fournie en préambule pour minimiser $f$ en partant de x0=(0,...,0). Tracer sur une courbe l'évolution de la valeur de la fonction au cours des itérations.

In [ ]:
# à compléter

Dans la suite, on mets en oeuvre un algorithme de descente coordonnée par coordonnée.

QN5 Écrire une fonction $T(x,i)$ calculant le minimum du problème (5) $$ \min_{t\in \mathbb{R}} f(x + t e_i),$$ où $e_i$ est le $i$ème vecteur de la base canonique. Pour ce faire,

  • on utilisera la question 6 de la feuille de TP pour réécrire $g:t\mapsto f(x+t e_i)$ sous la forme $g(t) = a t^2 + b t + \sqrt{\varepsilon + (c+t)^2} +
      D$
  • la fonction python newton_pure pour calculer le minimiseur de $g$.
In [ ]:
# à compléter

L'algorithme de descente coordonnée par coordonnée s'écrit sous la forme suivante. Pour $k\in \mathbb{N}$, on notera $i_k = 1 + \mathrm{mod}(k,d) \in \{1,\dots,d\}$, où $\mathrm{mod}(k,d)$ est le reste de la division euclidienne de $k$ par $d$. On considère les itérées $x^{(k)}$ définies par

$$ x^{(0)}\in\mathbb{R}^d,\quad x^{(k+1)} = x^{(k)} + T(x^{(k)}, i_k) e_{i_k} $$

QN6 Tracer $f(x^{(k)}) - $ fmin en échelle logarithmique, où $1\leq k\leq 2000$ et la suite $x^{(k)}$ est définie par la formule ci-dessus. On rappelle que pour calculer le reste de la division euclidienne de $k$ par $d$ on peut écrire k % d.

In [ ]:
# à compléter