Algorithmes d'optimisation -- L3 MINT et doubles licences -- Université Paris-Sud


$\newcommand{\Rsp}{\mathbb{R}} \newcommand{\nr}[1]{\|#1\|} \newcommand{\abs}[1]{|#1|} \newcommand{\eps}{\varepsilon} \newcommand{\sca}[2]{\langle#1|#2\rangle} \newcommand{\D}{\mathrm{D}} \newcommand{\hdots}{\dots} \newcommand{\cond}{\mathrm{cond}}$

TP 2: Régression logistique

Partie I: Descente de gradient avec rebroussement d'Armijo

Étant donné un point $x^{(k)}$ et une direction de descente $d^{(k)}$ le pas d'Armijo est défini de la manière suivante:

$$ t^{(k)}_{\mathrm{arm}} = \max \left\{ t \mid t = \beta^\ell, \ell\in \mathbb{N}, \hbox{ tq } f(x^{(k)} + t d^{(k)}) \leq f(x) + \alpha t\sca{\nabla f(x)}{d^{(k)}} \right\}, $$

où on prendra $\alpha = 0.3, \beta = 0.5$. Algorithmiquement, on procède de la manière suivante pour calculer $t^{(k)}_{\mathrm{arm}}$ en fonction de $x=x^{(k)}$ et $v = d^{(k)}$:

$$ \left| \begin{aligned} &\mathrm{pas\_armijo}(f,x,v)\\ &\qquad t \gets 1\\ &\qquad m \gets \sca{\nabla f(x)}{v}\\ &\qquad \textbf{while } f(x + t v) > f(x) + \alpha t m \\ &\qquad\qquad t \gets \beta t\\ &\qquad \mathrm{return}~t \end{aligned}\right. $$

L'algorithme de descente de gradient avec rebroussement d'Armijo s'écrit alors de la manière suivante:

$$ \begin{cases} d^{(k)} = - \nabla f(x^{(k)}) \\ t^{(k)} = \mathrm{pas\_armijo}(f,x^{(k)}, d^{(k)})\\ x^{(k+1)} = x^{(k)} + t^{(k)} d^{(k)} \end{cases}$$

QI.1) Écrire les fonctions suivantes:

  • pas_armijo(f,gradf,x) prend en entrée $f:\Rsp^N\to\Rsp$, gradf$=\nabla f:\Rsp^N\to\Rsp^N$, $x,d\in\Rsp^N$. On choisira alpha$=0.3$, beta$=0.5$ et on retournera le pas $t$.
  • gradient_armijo(f,gradf,x0,err=1e-6,maxiter=1000) prenant en argument la fonction f, gradf, x0=$x^{(0)}$ et un critère d'arrêt err$\in\Rsp$, et qui effectue l'algorithme de descente de gradient avec pas d'Armijo en arrêtant la boucle dès que $\nr{d^{(k)}} \leq$ err où dès que $k\geq$ maxiter. Cette fonction retournera le point $x$ trouvé, un vecteur contenant la suite $f(x^{(k)})$ et un second vecteur contenant $\nr{\nabla f(x^{(k)})}$. Indications:
    • Pour la fonction gradient_armijo, ne pas hésiter à copier-coller gradient_optimal et à faire les modifications nécessaires
    • Ne stocker que l'itération courante et pas la suite des itérées
In [1]:
# on importe les modules numpy et pyplot
import numpy as np
import matplotlib.pyplot as plt
# les deux commandes suivante paramètrent l'affichage des figures
%matplotlib inline
plt.rcParams['figure.figsize'] = [9.,6.]

# les paramètres alpha et beta comportent une valeur par défaut, c'est-à-dire que 
# pas_armijo(f,x,d,m) est équivalent à pas_armijo(f,x,d,m,0.3,0.5)
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
    
# la descente de gradient prend en argument:
# f = la fonction à évaluer
# gradf = le gradient de f
# x0 = le point de départ
def gradient_armijo(f,gradf,x0,err=1e-6,maxiter=1000):
    # ne pas hésiter à copier-coller et adapter gradient_optimal...
    x = x0
    E = []
    F = []
    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, 
        d = -gradf(x)
        F.append(f(x))
        E.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(F),np.array(E)

QI.2) Tester la fonction gradient_armijo en dimension $N=1$, $f(x) = x^{2} + e^{x}$, $x_0 = 1$ et en prenant pour critère d'arrêt err$= 10^{-8}$. Vous devriez trouver $x^* = -0.3517...$.

In [2]:
# test de la fonction pas_armijo
f = lambda x: x**2 + np.exp(x)
gradf = lambda x: 2*x + np.exp(x)

x,E,F = gradient_armijo(f,gradf,x0=1,err=1e-8)
print('solution: %g\nnombre d\'itérations: %d' % (x, len(E)))
solution: -0.351734
nombre d'itérations: 19

Partie II: Régression logistique avec rebroussement d'Armijo en dimension $N=2$

L'explication de la méthode est à lire dans la feuille de TD. On commence par construire les données $(x_a)_{a \in A} \subseteq \Rsp^2$ où :

  • $A = \{1,\hdots,n\} = A_0\sqcup A_1$, $A_0 = \{0, \hdots, m-1,\}$, $A_1 = \{m, 2m-1\}$, $n = 2m$
  • Les points $(x_a)_{a\in A_0}$ sont tirés selon une gaussienne centrée en $(0.5,0.5)$
  • Les points $(x_a)_{a\in A_1}$ sont tirés selon une gaussienne centrée en $(-0.5,-0.5)$.

On cherche alors à construire, par optimisation de la fonction $F$, un vecteur $w\in\Rsp^2$ tel que $\sigma(\sca{w}{x_a}) \simeq 0$ si $a\in A_0$ et $\sigma(\sca{w}{x_a})\simeq 1$ si $a\in A_1$.

In [3]:
m = 30
n = 2*m
X = np.vstack((np.hstack((-0.5+.2*np.random.randn(m,1), -0.5 + .2*np.random.randn(m,1))),
               np.hstack((0.5+.2*np.random.randn(m,1), 0.5 + .2*np.random.randn(m,1)))))
A0 = range(0,m)
A1 = range(m,2*m)
plt.plot(X[A0,0],X[A0,1],'.r',label='A0')
plt.plot(X[A1,0],X[A1,1],'.g',label='A1')
plt.legend()
plt.axis('equal');

QII.1: Écrire deux fonctions F(w,X,A0,A1,gamma) et gradF(w,X,A0,A1,\gamma) calculant la valeur et le gradient de la fonctionnelle

$$ F(w) = - \left(\sum_{a\in A_0} \log(1-\sigma(\sca{w}{x_a})) + \sum_{a\in A_1} \log(\sigma(\sca{w}{x_a})\right) + \frac{\gamma}{2} \nr{w}^2. $$

Tester que gradF correspond bien au gradient de F en utilisant la fonction verifier_gradient fournie ci-dessous.

In [4]:
# on définit une fonction sigma(t) correspondant à la sigmoïde
def sigma(t):
    return 1/(1+np.exp(-t))

def F(w,X,A0,A1,gamma):
    # on calcule F en utilisant la formule (on pourra utiliser np.sum pour gagner en concision/rapidité)
    h = sigma(np.dot(X,w.T))
    return -(np.sum(np.log(1.-h[A0])) + np.sum(np.log(h[A1]))) + .5*gamma*np.dot(w,w)

def gradF(w,X,A0,A1,gamma):
    # on calcule le gradient de F en appliquant la formule obtenue en TD
    d = len(w)
    g = np.zeros(d)
    h = sigma(np.dot(X,w.T))
    for i in A0: g += h[i]*X[i,:]
    for j in A1: g += (h[j]-1)*X[j,:]
    g += gamma*w
    return g

""" la fonction verifier_gradient calcule le gradient de f en un point x0 de deux manières:
    1) via la fonction gradf donnée en argument 
    2) via des différences finies: df/de_i (x_0) ~= (f(x_0+eps e_i) - f(x_0))/eps
    et compare ces deux résultats: si la fonction gradf est correctement programmée, la différence
    entre les deux méthode de calcul devrait être petite """
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))
    
# on vérifie l'implémentation de gradF en utilisant la fonction verifier_gradient
gamma = .1
f = lambda w: F(w,X,A0,A1,gamma)
gradf = lambda w: gradF(w,X,A0,A1,gamma)
verifier_gradient(f, gradf, np.random.rand(2))
erreur numérique dans le calcul du gradient: 2.11816e-10 (doit être petit)

QI.2: Calculer $w=\arg\min_{\Rsp^2} F$ en utilisant la fonction gradient_armijo. Tracer sur un graphe l'évolution de la norme du gradient au cours des itérations de la descente de gradient, soit $(k,\nr{\nabla F(w^{(k)})})$, où l'axe $y$ sera en échelle logarithmique.

In [5]:
w0 = np.zeros(2) 
gamma = .01
w,fiter,giter = gradient_armijo(f, gradf, w0, err=1e-7)
plt.semilogy(giter)
Out[5]:
[<matplotlib.lines.Line2D at 0x7ff260d36a58>]

QI.3: Soit $w$ le minimiseur approché de $F$ calculé dans la question précédente. Tracer sur un graphe:

  • les points $(x_a)_{a\in A_0}$ en rouge
  • les points $(x_a)_{a\in A_1}$ en bleu
  • quelques lignes de niveau de la fonction $u:x\mapsto \sigma(\sca{w}{x})$ en utilisant la fonction fournie afficher_lignes_niveau(w).

Faire varier le paramètre de régularisation $\gamma \in \{0.01,0.1,1,10\}$ et observer:

  • comment varie le nombre d'itérations.
  • comment varie $\nr{w}$ (où $\nr{w}$ est le minimiseur).
  • quel est l'effet sur la fonction $u:x\mapsto \sigma(\sca{w}{x_i})$.
In [6]:
# on définit une fonction permettant d'afficher les lignes de niveaux de u_w: R^2 -> R
# definie par u_w(x) = \sigma(<x|w>)
# NB: u peut être définie en utilisant lambda, par exemple:
#    u = lambda (x,y): ... où x,y sont les coordonnées du point
def afficher_lignes_niveau(w):
    Xcontour,Ycontour = np.meshgrid(np.linspace(-1., 1., 100),
                                    np.linspace(-1., 1., 100))
    u = lambda x,y: sigma(x*w[0] + y*w[1])
    Zcontour = u(Xcontour,Ycontour)
    plt.axes([-1.5, -1.5, 1.5, 1.5])
    plt.axis('equal')
    p = plt.contour(Xcontour, Ycontour, Zcontour, cmap='RdBu')
    plt.clabel(p, inline=1, fontsize=10)
    
for gamma in [0.01,0.1,1,5]:
    w,fiter,giter = gradient_armijo(lambda w: F(w,X,A0,A1,gamma),
                                    lambda w: gradF(w,X,A0,A1,gamma), 
                                    np.zeros(2), err=1e-7)
    print('gamma = %d: niter = %d, |w|=%g' % (gamma,len(fiter),np.linalg.norm(w)))
    plt.figure()
    afficher_lignes_niveau(w)
    plt.scatter(X[A0,0],X[A0,1],marker='s',color='r', label='(x_a), a\in A_0')
    plt.scatter(X[A1,0],X[A1,1],marker='o',color='b', label='(x_a), a\in A_1')
    plt.axis('tight')
    plt.axis('equal')
    plt.title('donn\\\'ees $(x_i)_{i\in I\cup J}$ et solution $u: x \\mapsto \sigma(\\langle x | w\\rangle)$ / ')
    plt.legend()
gamma = 0: niter = 495, |w|=9.8531
gamma = 0: niter = 74, |w|=6.24776
gamma = 1: niter = 16, |w|=3.39274
gamma = 5: niter = 17, |w|=1.78803

Partie III: Régression logistique avec méthode de Newton en dimension $N=2$

QIII.1: Écrire une fonction hessF calculant la hessienne de $F$, en utilisant la formule obtenue en TD. Vérifier le bon fonctionnement de hessF via la fonction verifier_hessienne fournie ci-dessous.

In [7]:
""" cette fonction vérifie numériquement que hessf (=fonction évaluant la hessienne de f) est bien
    la dérivée de grad (= fonction évaluant le gradient de f), en un point x0 """
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)
    
def hessF(w,X,A0,A1,K):
    h = sigma(np.dot(X,w.T))
    d = len(w) 
    N = X.shape[0]
    A = K*np.eye(d)
    for i in range(N):
        A += h[i] * (1-h[i]) * np.tensordot(X[i,:],X[i,:],axes=0)    
    return A

# on vérifie que la hessienne est correcte en utilisant la fonction verifier_hessienne
gamma = .1
f = lambda w: F(w,X,A0,A1,gamma)
gradf = lambda w: gradF(w,X,A0,A1,gamma)
hessf = lambda w: hessF(w,X,A0,A1,gamma)
verifier_hessienne(gradf, hessf, x0=np.random.rand(2))
erreur numerique dans le calcul de la hessienne: 1.92549e-19 (doit etre petit)

QIII.2: Écrire une fonction newton_armijo(f,gradf,hessf,x0,err=1e-6,maxiter=1000) implémentant la méthode de Newton avec rebroussement d'Armijo, c'est-à-dire l'algorithme suivant:

$$ \begin{cases} g^{(k)} = \nabla F(x^{(k)}) \\ d^{(k)} = -[D^2 F(x^{(k)})]^{-1} \nabla F(x^{(k)}) \\ t^{(k)} = \mathrm{pas\_armijo}(f,x^{(k)},d^{(k)})\\ x^{(k+1)} = x^{(k)} + t^{(k)} d^{(k)} \end{cases} $$

Il ne faut pas hésiter à copier-coller la fonction gradient_armijo et à changer ce qui est nécessaire. Pour la résolution du système linéaire, on utilisera np.linalg.solve.

Utiliser cette fonction pour calculer le minimum de $F$. Comme dans la question Q1.2, tracer la norme du gradient en fonction de l'itération. Vérifier également que le nombre d'itérations dépend peu du choix du paramètre $\gamma$.

In [8]:
def newton_armijo(f,gradf,hessf,x0,err=1e-6,maxiter=500):
    # la boucle est tres similaire à celle de gradient_armijo: la seule chose qui
    # change est le choix de la direction de descente. 
    # NB: on utilise np.linalg.solve pour résoudre le système linéaire apparaissant dans la méthode
    x = x0.copy()
    fIter = []
    gIter = []
    k = 0
    while(True):
        k = k+1
        if k > maxiter: # maximum de 10^6 itérations
            print('erreur: nombre maximum d\'itérations atteint')
            break
        gx = -gradf(x)
        dx = np.linalg.solve(hessf(x),gx)
        fIter.append(f(x))
        gIter.append(np.linalg.norm(dx))
        if np.abs(np.linalg.norm(gx)) <= err:
            break
        t = pas_armijo(f,gradf,x,dx)
        print('iteration %d: f=%f, |g|=%f, step=%f' % (k, f(x), np.linalg.norm(gx),t))
        x = x + t*dx
    return x,np.array(fIter),np.array(gIter)


# comparer le nombre d'iterations dans la methode de Newton pour gamma dans [0.01,0.1,1,5]
for gamma in [0.01,0.1,1,10]:
    w0 = np.zeros(2)
    print("gamma=%g"%gamma)
    w,fiter,giter = newton_armijo(lambda w: F(w,X,A0,A1,gamma),
                                  lambda w: gradF(w,X,A0,A1,gamma),
                                  lambda w: hessF(w,X,A0,A1,gamma), w0, 1e-7, 1000)
    plt.semilogy(giter,label=('gamma=%g'%gamma))
plt.legend()
gamma=0.01
iteration 1: f=41.588831, |g|=21.745132, step=1.000000
iteration 2: f=9.413018, |g|=5.530196, step=1.000000
iteration 3: f=3.732362, |g|=2.004350, step=1.000000
iteration 4: f=1.697714, |g|=0.736485, step=1.000000
iteration 5: f=0.964274, |g|=0.257311, step=1.000000
iteration 6: f=0.749175, |g|=0.074768, step=1.000000
iteration 7: f=0.714942, |g|=0.012777, step=1.000000
iteration 8: f=0.713598, |g|=0.000563, step=1.000000
iteration 9: f=0.713595, |g|=0.000001, step=1.000000
gamma=0.1
iteration 1: f=41.588831, |g|=21.745132, step=1.000000
iteration 2: f=9.861540, |g|=5.393788, step=1.000000
iteration 3: f=4.760083, |g|=1.809644, step=1.000000
iteration 4: f=3.398553, |g|=0.518892, step=1.000000
iteration 5: f=3.188201, |g|=0.086899, step=1.000000
iteration 6: f=3.180330, |g|=0.003597, step=1.000000
iteration 7: f=3.180316, |g|=0.000007, step=1.000000
gamma=1
iteration 1: f=41.588831, |g|=21.745132, step=1.000000
iteration 2: f=13.705947, |g|=4.239893, step=1.000000
iteration 3: f=11.674223, |g|=0.711389, step=1.000000
iteration 4: f=11.597220, |g|=0.030901, step=1.000000
iteration 5: f=11.597065, |g|=0.000064, step=1.000000
gamma=10
iteration 1: f=41.588831, |g|=21.745132, step=1.000000
iteration 2: f=28.519309, |g|=0.719227, step=1.000000
iteration 3: f=28.503780, |g|=0.002261, step=1.000000
Out[8]:
<matplotlib.legend.Legend at 0x7ff260c9bd30>

Partie IV: Application à la classification de caractères manuscrits

L'objectif de cette partie est d'appliquer la méthode décrite dans les parties précédentes (où l'on traitait des données en dimension $N=2$) à un jeu de données plus réaliste et de dimension beaucoup plus grande. L'objectif est de distinguer des images de caractères manuscrits. Pour cela, on a demandé à 6000 personnes d'écrire les chiffres de zéro à neuf. Ici, l'on cherchera uniqument à distinguer les chiffres zéro manuscrits des chiffres un. Chaque point de donnée $x_a$ est construit de la manière suivante:

  • On demande à une personne d'écrire un chiffre (zéro, un, ..., neuf), que l'on stocke sous la forme d'une image en niveaux de gris de taille 28x28 pixels. On représente cette image par une matrice carrée $p_a \in \mathcal{M}_{28,28}(\Rsp)$ (où $(p_{a})_{ij} = 0$ si le pixel $(i,j)$ est noir et $1$ s'il est blanc).
  • La matrice carrée $p_a \in \mathcal{M}_{28,28}(\Rsp)$ est convertie en un vecteur $x_a\in\Rsp^{N}$ de taille $N = 28\times 28 = 784$. Pour cela, on met bout-à-bout les lignes de la matrice $p$, via la fonction np.reshape.
  • Enfin, on stocke séparément un entier correspondant au chiffre manuscrit (on appelle ça le "label", ou étiquette, de l'image).

Ce faisant, on obtient donc deux tableaux:

  • images possède $60000$ lignes (soit $10$ chiffres zéro, un, $\hdots$ dessinés par $6000$ personnes) et $784=28^2$ colonnes (images de $28\times 28$ pixels stockées sous formes de vecteurs). images[i] correspond alors à la $i$ème image, représentée par un vecteur de $\Rsp^{768}$.
  • labels possède $60000$ lignes et une seule colonne: label[i] contient le numéro du chiffre manuscrit représenté par l'image images[i]

La fonction read_mnist() ci-dessous télécharge ces données et retourne les tableaux images et labels. On fournit une fonction d'affichage: show_mnist(images[i]) affiche la $i$ème image.

In [14]:
# import du jeu de données MNIST de Yan Lecun: 60000 chiffres manuscrits avec labels
import os
import struct
import urllib
import gzip
import shutil
try:
    from urllib.request import urlretrieve
except ImportError:
    from urllib import urlretrieve

def gunzip(iname,oname):
    with gzip.open(iname, 'rb') as f_in:
        with open(oname, 'wb') as f_out:
            shutil.copyfileobj(f_in, f_out)
def download_mnist(dest):
    src = 'http://yann.lecun.com/exdb/mnist'
    for f in ['train-images-idx3-ubyte', 'train-labels-idx1-ubyte']:
        urlretrieve('%s/%s.gz' % (src,f), '%s/%s.gz' % (dest,f))
        gunzip('%s/%s.gz' % (dest,f),'%s/%s' % (dest,f))
def read_mnist():
    d = os.path.expanduser('~/.m315-mnist-data')
    if not os.path.isdir(d):
        os.makedirs(d)
        download_mnist(d)
    with open('%s/train-labels-idx1-ubyte' % d, 'rb') as flbl:
        magic, num = struct.unpack(">II", flbl.read(8))
        labels = np.fromfile(flbl, dtype=np.int8)
    with open('%s/train-images-idx3-ubyte' % d, 'rb') as fimg:
        magic, num, rows, cols = struct.unpack(">IIII", fimg.read(16))
        images = np.fromfile(fimg, dtype=np.uint8).reshape(len(labels), rows*cols)/255.
    return images,labels
def show_mnist(image, cmap='gray'):
    plt.figure()
    imgplot = plt.imshow(np.reshape(image,(28,-1)), interpolation='nearest', cmap=cmap)
   
images,labels = read_mnist()

QIV.0 Passer quelques minutes à jouer avec les fonctions show_mnist() et print() (afficher les labels, et aussi les vecteurs images[i]) pour comprendre le jeu de données.

In [15]:
# exemple d'image: faire varier un peu i pour comprendre le jeu de données
i = 16
x = images[i]
p = np.reshape(images[i],(28,28))
show_mnist(p)
print(x.shape)
print(p.shape)
print(labels[i])
(784,)
(28, 28)
2

QIV.1 Construire un tableau X0 contenant $m=100$ images du chiffre zéro extraite du tableau images, et un tableau X1 contenant également $m$ représentants du chiffre un également extraits de images (X0 et X1 possèderont donc $m$ lignes et $784$ colonnes). Construire ensuite les données du problème d'optimisation:

  • Un tableau X de $n=2m$ lignes et $784$ colonnes obtenu en mettant à bout X0 et X1 (on pourra utiliser np.vstack).
  • Un tableau A0 (resp. A1) contenant la liste des indices correspondant aux chiffres zéro (resp un).
In [16]:
ntot = len(labels)
m = 100 # nombre d'images par classe
X0 = [images[i] for i in range(ntot) if labels[i] == 0][0:m]
X1 = [images[i] for i in range(ntot) if labels[i] == 1][0:m]
X = np.vstack((X0,X1))
A0 = range(0,m)
A1 = range(m,2*m)

gamma = .1
f = lambda w: F(w,X,A0,A1,gamma)
gradf = lambda w: gradF(w,X,A0,A1,gamma)
hessf = lambda w: hessF(w,X,A0,A1,gamma)

verifier_gradient(f, gradf, x0=np.zeros(784))
verifier_hessienne(gradf, hessf, x0=np.random.rand(784))
erreur numérique dans le calcul du gradient: 2.14485e-08 (doit être petit)
erreur numerique dans le calcul de la hessienne: 6.43181e-14 (doit etre petit)

QIV.2 Résoudre le problème d'optimisation $\min_{w\in\Rsp^d} F(w)$ pour $\gamma=0.1$ en utilisant la méthode de Newton.

In [12]:
gamma = .1
w0 = np.zeros(784)
w,fIter,gIter = newton_armijo(lambda w: F(w,X,A0,A1,gamma),
                              lambda w: gradF(w,X,A0,A1,gamma),
                              lambda w: hessF(w,X,A0,A1,gamma),w0,
                              err=1e-8,maxiter=5000)
iteration 1: f=138.629436, |g|=377.194291, step=1.000000
iteration 2: f=25.816351, |g|=89.958706, step=1.000000
iteration 3: f=9.159966, |g|=31.569106, step=1.000000
iteration 4: f=3.712922, |g|=11.593728, step=1.000000
iteration 5: f=1.701434, |g|=4.355502, step=1.000000
iteration 6: f=0.926405, |g|=1.640519, step=1.000000
iteration 7: f=0.640555, |g|=0.598944, step=1.000000
iteration 8: f=0.557049, |g|=0.190618, step=1.000000
iteration 9: f=0.543141, |g|=0.041410, step=1.000000
iteration 10: f=0.542359, |g|=0.003433, step=1.000000
iteration 11: f=0.542354, |g|=0.000029, step=1.000000

Étant donné $w^* \in \Rsp^{768}$ l'unique solution du problème d'optimisation, on notera $u^*: x\mapsto\sigma(\sca{x}{w^*})$. Le principe de la régression logistique (cf intro du TD) est que $u^*(x_a) \simeq 1$ si $a \in A_0$ (c'est-à-dire si l'image $x_a$ représente le chiffre zéro) et $u^*(x_a) \simeq 0$ si $a\in A_1$ (si l'image $x_a$ représente le chiffre un). Comme on a utilisé seulement une petite fraction des images pour construire $w$ (on a choisi $n=100$ images du chiffre zéro alors que la base de données en contient $6000$!), on va maintenant chercher à valider la méthode en comparant la catégorie ($0$ où $1$) devinée par la fonction $u^*$ à celle enregistrée dans le tableau labels sur l'ensemble des zéros et des uns manuscrits contenus dans la base de données.

QIV.3 Vérifier la pertinence de la méthode de la manière suivante :

  • Pour chaque image MNIST $x = images[i]$ représentant un $0$ où un $1$ (pour tout $i$ tel que labels[i] $\in \{0,1\}$), calculer $u^*(x)$, et définir label_devine=0 si $u(x)< 1/2$ et label_devine=$1$ sinon.
  • Comparer label_devine au véritables label (i.e. labels[i]).
  • Calculer enfin la proportion de cas où label_devine=labels[i]: ce nombre mesure la proportion de cas où la régression logistique permet de reconnaître correctement le chiffre représenté par l'image.
In [108]:
total = 0
correct = 0
for i in range(ntot):
    if labels[i] not in [0,1]:
        continue
    P = sigma(np.dot(images[i],w.T)) > .5
    total = total+1
    if P==labels[i]:
        correct = correct+1
print('pourcentage de cas ou le chiffre 0 est correctement reconnu: %g%%' % (100.*float(correct)/float(total)))
pourcentage de cas ou le chiffre 0 est correctement reconnu: 95.0336%

QIV.4: Afficher le vecteur $w$ sous la forme d'une image en utilisant show_mnist(). Interpréter.

In [104]:
show_mnist(np.reshape(w,(28,28)),cmap='RdBu')

QIV.5 Adapter l'approche à la séparation des chiffres 0 de tout les autres chiffres. (on pourra prendre m=100 chiffres $0$ et $m$ chiffres de chaque autre classe, soit 1000 points de donnée au total).

In [107]:
n = 100 # nombre d'images par classe
X0 = [images[i] for i in range(ntot) if labels[i] == 0][0:n]
X1 = [images[i] for i in range(ntot) if labels[i] >= 1][0:9*n]
X = np.vstack((X0,X1))
A0 = range(0,n)
A1 = range(n,10*n)

gamma = 10
w0 = np.zeros(784)
w,fIter,gIter = newton_armijo(lambda w: F(w,X,A0,A1,gamma),
                              lambda w: gradF(w,X,A0,A1,gamma),
                              lambda w: hessF(w,X,A0,A1,gamma),w0, 
                                    err=1e-8,maxiter=5000)
show_mnist(np.reshape(w,(28,28)),cmap='RdBu')

total = 0
correct = 0
for i in range(ntot):
    R = labels[i] >= 1
    P = sigma(np.dot(images[i],w.T)) > .5
    total = total+1
    if P==R:
        correct = correct+1
print('pourcentage de cas ou le chiffre 0 est correctement reconnu: %g%%' % (100.*float(correct)/float(total)))
iteration 1: f=693.147181, |g|=2277.565327, step=1.000000
iteration 2: f=186.520113, |g|=572.343704, step=1.000000
iteration 3: f=98.207586, |g|=205.937616, step=1.000000
iteration 4: f=69.212164, |g|=71.872824, step=1.000000
iteration 5: f=61.692959, |g|=20.282336, step=1.000000
iteration 6: f=60.731338, |g|=3.026784, step=1.000000
iteration 7: f=60.706153, |g|=0.092316, step=1.000000
iteration 8: f=60.706130, |g|=0.000092, step=1.000000
pourcentage de cas ou le chiffre 0 est correctement reconnu: 98.6%