** Algorithmes d'optimisation -- L3 MINT et doubles licences 2017/2018 -- 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}}$
L'explication de la méthode est à lire dans la feuille de TD. On commence par construire les données $(x_i)_{i\in I \cup J} \subseteq \Rsp^2$:
On cherche alors à construire, par optimisation de la fonction $F$, un vecteur $w\in\Rsp^2$ tel que $\sigma(\sca{w}{x_i}) \simeq 0$ si $i\in I$ et $\sigma(\sca{w}{x_i})\simeq 1$ si $i\in J$.
On définit également des fonctions tirées des TP précédents:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# la commande suivante agrandit les figures
plt.rcParams['figure.figsize'] = [9.,6.]
# initialiser les variables K,Q, x0, fmin
n = 30
X = np.vstack((np.hstack((-0.5+.2*np.random.randn(n,1), -0.5 + .2*np.random.randn(n,1))),
np.hstack((0.5+.2*np.random.randn(n,1), 0.5 + .2*np.random.randn(n,1)))))
I = range(0,n)
J = range(n,2*n)
plt.plot(X[I,0],X[I,1],'.r')
plt.plot(X[J,0],X[J,1],'.g')
plt.axis('equal');
# on recupere ici des fonctions issues des TP précédents. Seule la fonction check_hessian est nouvelle.
def backtrack(f,x,d,m,alpha=0.3,beta=0.5):
t = 1
while f(x+t*d) > f(x) + alpha*t*m:
t = beta*t
return t
def gradient_backtracking(f,g,x0,err=1e-6,maxiter=500):
x = x0.copy()
fiter = []
giter = []
k = 0 # nombre d'itérations
while(True):
k = k+1
if k > maxiter: # maximum de 10^6 itérations
print('erreur: nombre maximum d\'itérations atteint')
break
d = -g(x)
fiter.append(f(x))
giter.append(np.linalg.norm(d))
if np.linalg.norm(d) <= err:
break
t = backtrack(f,x,d,-np.linalg.norm(d)**2)
if k%10==0: # on affiche des informations toute les 20 itérations
print('iteration %d: f=%f, |g|=%f, step=%f' % (k, f(x), np.linalg.norm(d),t))
x = x + t*d
return x,np.array(fiter),np.array(giter)
def check_gradient(f,g,x0):
N = len(x0)
gg = np.zeros(N)
for i in range(N):
eps = 1e-4
e = np.zeros(N)
e[i] = eps
gg[i] = (f(x0+e) - f(x0-e))/(2*eps)
print('erreur numerique dans le calcul du gradient: %g (doit etre petit)' % np.linalg.norm(g(x0)-gg))
# cette fonction vérifie que h (=fonction évaluant le hessian) est bien la dérivée
# de g (= fonction évaluant le gradient), en un point x0
def check_hessian(g,h,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,:] = (g(x0+e) - g(x0-e))/(2*eps)
print('erreur numerique dans le calcul de la hessienne: %g (doit etre petit)' % np.sum((H-h(x0)))**2)
QI.1: Écrire une fonction $F(w,X,I,J,\gamma)$ et $gF(w,X,I,J,\gamma)$ calculant la valeur et le gradient de la fonctionnelle
$$ F(w) = - \left(\sum_{i\in I} \log(1-\sigma(\sca{w}{x_i})) + \sum_{i\in J} \log(\sigma(\sca{w}{x_i})\right) + \frac{\gamma}{2} \nr{w}^2. $$Tester que $gF$ correspond bien au gradient de $F$ en utilisant la fonction check_gradient.
# on définit une fonction sigma(t) correspondant à la sigmoïde
# <completer>
def F(w,X,I,J,gamma):
# on calcule F en utilisant la formule (on pourra utiliser np.sum)
# <completer>
def gF(w,X,I,J,gamma):
d = len(w)
g = np.zeros(d)
# on calcule le gradient en appliquant la formule obtenue en TD
# <completer>
return g
# on vérifie le calcul du gradient en utilisant la fonction check_gradient
# rappel: on peut utiliser f = lambda w: F(w,X,I,J,gamma) pour définir
# f(w) = F(w,X,I,J,gamma)
# <completer>
QI.2: Calculer $w=\arg\min_{\Rsp^2} F$ en utilisant la fonction gradient_backtracking. 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.
# <completer>
QI.3: Soit $w$ le minimiseur approché de $F$ calculé dans la question précédente. Tracer sur un graphe:
Faire varier le paramètre de régularisation $\gamma \in \{0.01,0.1,1,10\}$ et observer:
# on définit une fonction permettant d'afficher les lignes de niveaux de u: R^2 -> R
# 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 lignes_niveau(u):
Xcontour,Ycontour = np.meshgrid(np.linspace(-1., 1., 100),
np.linspace(-1., 1., 100))
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)
# <completer>
QI.4: Écrire une fonction $hF$ calculant la hessienne de $F$, en utilisant la formule du TD. Vérifier son bon fonctionnement via check_hessian (voir l'introduction du TP).
# <completer>
Q1.5: Écrire une fonction newton_backtracking implémentant la méthode de Newton avec backtracking (ne pas hésiter à copier-coller la fonction gradient_backtracking et à changer ce qui est nécessaire) :
$$ \begin{cases} g^{(k)} = \nabla F(x^{(k)}) \\ d^{(k)} = -D^2 F(x^{(k)})^{-1} \nabla F(x^{(k)}) \\ m = \sca{d^{(k)}}{g^{(k)}} \\ t^{(k)} = \hbox{obtenu par backtracking d'Armijo, cf la fonction backtrack}\\ x^{(k+1)} = x^{(k)} + t^{(k)} d^{(k)} \end{cases} $$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$.
def newton_backtracking(f,g,h,x0,err=1e-6,maxiter=500):
# la boucle est tres similaire à celle de gradient_backtracking: 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
# <completer>
# comparer le nombre d'iterations dans la methode de Newton pour gamma dans [0.01,0.1,1,10]
# <completer>
On applique la méthode décrite dans la partie I (où l'on traitait des données en dimension $d=2$) à un jeu de données de dimension beaucoup plus grande. L'objectif est de distinguer des caractères manuscrits: pour cela, on a demandé à 6000 personnes d'écrire les chiffres de $0$ à $9$. Ici, l'on cherchera à distinguer les $0$ des $1$. Chaque point $x_i \in \Rsp^{784}$ est construit de la manière suivante:
La fonction read_mnist() ci-dessous télécharge ces données et retourne deux tableaux numpy images et labels:
# MNIST
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,28)), interpolation='nearest', cmap=cmap)
images,labels = read_mnist()
Il faut 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.
# exemple d'image: faire varier un peu i pour comprendre le jeu de données
i = 15
show_mnist(images[i])
print(labels[i])
7
QII.1 Construire un tableau $X0$ contenant $n=100$ images du chiffre $0$ extraite du tableau 'images', et un tableau $X1$ contenant également $n$ représentants du chiffre $1$ également extraits de images (X0 et X1 possèderont donc $n$ lignes et $784$ colonnes). Construire ensuite les données du problème d'optimisation:
# <completer>
QII.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.
# <completer>
Soit $w$ la solution du problème d'optimisation et $u:x\mapsto\sigma(\sca{x}{w})$. L'idée de la régression logistique (cf intro du TD) est que $u(x) \simeq 1$ si $x$ représente un $1$ et $u(x) \simeq 0$ si $x$ représente un $0$. On n'a utilisé qu'une petite fraction des images pour construire $w$ (on a choisi $n=100$ images de $0$ alors que la base de données en contient $6000$): l'idée est de valider la méthode en comparant la classe ($0$ où $1$) devinée par la fonction $u$ à celle enregistrée dans le tableau labels sur l'ensemble des $0$ et $1$ de la base de données.
QII.3 Vérifier la pertinence de la méthode de la manière suivante :
total = 0
correct = 0
# <completer>
print('pourcentage de cas ou le chiffre 0 est correctement reconnu: %g%%' % (100.*float(correct)/float(total)))
QII.4: Afficher le vecteur $w$ sous la forme d'une image en utilisant show_mnist(). Interpréter.
# <completer>
QII.5 Adapter l'approche à la séparation des chiffres 0 de tout les autres chiffres. (on pourra prendre n=100 chiffres $0$ et $n$ chiffres de chaque autre classe, soit 1000 points de donnée au total).
# <completer>