In [2]:
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 [3]:
def g(a,b,c,eps,t):
    return a*t**2 + b*t + np.sqrt(eps + (c+t)**2)
def gprime(a,b,c,eps,t):
    return 2*a*t + b + (c+t)/np.sqrt(eps + (c+t)**2)
def gseconde(a,b,c,eps,t):
    return 2*a + eps/(eps + (c+t)**2)**1.5

## 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))
erreur numérique dans le calcul du gradient: 7.1334e-12 (doit être petit)
erreur numérique dans le calcul du gradient: 7.7543e-12 (doit être petit)

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 [4]:
def newton_pur(a, b, c, eps):
    t = -c
    for k in range(5):
        t = t - gprime(a,b,c,eps,t)/gseconde(a,b,c,eps,t)
    return t

## 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))
cas 1: solution trouvée -5.00039, solution attendue -5.0003905411
cas 2: solution trouvée -995, solution attendue -995.0000000024029

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 [5]:
# 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 [7]:
def f(x):
    return .5*mu*np.linalg.norm(A@x - y)**2 + np.sum(np.sqrt(eps + x**2))

def gradf(x):
    return mu*(A.T@[email protected] - A.T@y) + x/np.sqrt(eps+x**2)

## vérification du gradient de f
verifier_derivee(f,gradf,np.random.rand(d))
erreur numérique dans le calcul du gradient: 5.76011e-08 (doit être petit)

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 [8]:
x,F = gradient_armijo(f,gradf,x0=np.zeros(d))
plt.semilogy(F - fmin)
Out[8]:
[<matplotlib.lines.Line2D at 0x7f3c0c0aaef0>]

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 [9]:
def T(x,i):
    ei = np.zeros_like(x)
    ei[i] = 1
    a = .5*mu*np.linalg.norm(A@ei)**2
    b = mu*(A@x-y)@(A@ei)
    c = x[i]
    t = newton_pur(a, b, c, eps)
    return t

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 [10]:
niter = 2000
F = np.zeros(niter)
x = np.zeros(d)
for k in range(niter):
    i = k % d
    ei = np.zeros_like(x)
    ei[i] = 1
    x[i] += T(x,i)
    F[k] = f(x)
plt.semilogy(F - fmin)
Out[10]:
[<matplotlib.lines.Line2D at 0x7f3bf076ad30>]