Présentation du problème de reconnaissance de caractères manuscrits (MNIST DataBase à partir d’images numérisées. L’objectif est de comparer les performances (qualité de prévision, temps d'exécution) en fonction de latechnologie, ici Python et la librairie Scikit-learn, et en fonction de la taille de l'échantillon. Même interprété, les exécutions sont effecaces par une bonne parallélisation. La faiblesse concerne les insufisances des aides à l'interprétation..
L'objectif général est la construction d'un meilleur modèle de reconnaissance de chiffres manuscrits. Ce problème est ancien (zipcodes) et sert souvent de base pour la comparaison de méthodes et d'algorithmes d'apprentissage. Le site de Yann Le Cun: MNIST DataBase, est à la source des données étudiées, il décrit précisément le problème et les modes d'acquisition. Il tenait à jour la liste des publications proposant des solutions avec la qualité de prévision obtenue. Ce problème a également été proposé comme sujet d'un concours Kaggle mais sur un sous-ensemble des données.
De façon très schématique, plusieurs stratégies sont développées dans une vaste littérature sur ces données.
L'objectif de cet atelier est de comparer sur des données relativement volumineuses les performances de différents environnements technologiques et librairies. Une dernière question est abordée, elle concerne l'influence de la taille de l'échantillon d'apprentissage sur le temps d'exécution ainsi que sur la qualité des prévisions.
Analyse des données avec Python noter les temps d'exécution, la précision estimée sur l'échantillon test.
Les données peuvent être préalablement téléchargées ou directement lues. Ce sont celles originales du site MNIST DataBase mais préalablement converties au format .csv, certes plus volumineux mais plus facile à lire. Attention le fichier mnist_train.zip présent dans le dépôt est compressé.
In [ ]:
# Graphiques dans la fenêtre
%matplotlib inline
In [ ]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
Les données peuvent être préalablement téléchargées ou directement lues. Ce sont celles originales du site MNIST DataBase mais préalablement converties au format .csv, certes plus volumineux mais plus facile à lire. Attention le fichier mnist_train.zip présent dans le dépôt est compressé.
In [ ]:
# Lecture des données d'apprentissage
# path="" # Si les données sont dans le répertoire courant sinon:
path="http://www.math.univ-toulouse.fr/~besse/Wikistat/data/"
Dtrain=pd.read_csv(path+"mnist_train.csv",header=None)
Dtrain.head()
In [ ]:
# Extraction puis suppression de la dernière colonne des labels
Ltrain=Dtrain.iloc[:,784]
Dtrain.drop(Dtrain.columns[[784]], axis=1,inplace=True)
Dtrain.head()
In [ ]:
# Dimensions de l'échantillon
Dtrain.shape
In [ ]:
# Même chose pour les données de test
Dtest=pd.read_csv(path+"mnist_test.csv",header=None)
Ltest=Dtest.iloc[:,784]
Dtest.drop(Dtest.columns[[784]], axis=1,inplace=True)
Dtest.shape
In [ ]:
# affichage d'un chiffre
plt.figure(1, figsize=(3, 3))
plt.imshow(np.matrix(Dtest.iloc[1,:]).reshape(28,28), cmap=plt.cm.gray_r, interpolation='nearest')
plt.show()
Les données ont déjà été normalisées centrées et sont complètes. Elles ne nécessitent pas d'autre "nettoyage" au moins rudimentaire.
Le tutoriel d'introduction à Scikit-learn montre comment représenter les images des caractères ainsi qu'une ACP qui n'est pas reprise ici. Quelles sont néanmoins les performances de k-means sur un tel volume ?
In [ ]:
from sklearn.metrics import confusion_matrix
from sklearn.cluster import KMeans
tps1 = time.clock()
km=KMeans(n_clusters=10,init='k-means++',
n_init=10, max_iter=100, tol=0.01,
precompute_distances=True, verbose=0,
random_state=None, copy_x=True, n_jobs=1)
km.fit(Dtrain)
tps2 = time.clock()
print("Temps execution Kmeans :", (tps2 - tps1)/60)
In [ ]:
cm = confusion_matrix(Ltrain, km.labels_)
print(cm)
Résultat sans grand intérêt mais qui montre la difficulté de regouper les caractères identiques à l'aide de la distance euclidienne usuelle; il y a beaucoup de confusion entre les classes.
In [ ]:
# Définition du modèle avec un nombre k "standard" de voisins
from sklearn.neighbors import KNeighborsClassifier
tps1 = time.clock()
knn = KNeighborsClassifier(n_neighbors=10,n_jobs=-1)
digit_knn=knn.fit(Dtrain, Ltrain)
tps2 = time.clock()
print("Temps de k-nn :",(tps2 - tps1)/60)
In [ ]:
# Apprentissage et estimation de l'erreur de prévision sur l'échantillon test
tps1 = time.clock()
erreur=1-digit_knn.score(Dtest,Ltest)
tps2 = time.clock()
print("Temps:",(tps2 - tps1)/60,"Erreur:",erreur)
Il faudrait ré-appliquer la procédure d'otpimisation de $k$ par validation croisée décrite dans le tutoriel d'introduction à scikit-learn. Néanmoins la solution $k=10$ est raisonnable et on retrouve une performance classique sur ce type de données: 3.3%, pour une méthode utilisée sans raffinement.
C'est en effet une autre distance qu'il faudrait utiliser avec les $k$ plus proches voisins pour améliorer sensiblement les résultats mais avec un coût beaucoup plus élevé en temps de calcul. Un autre scénario propose ainsi le calcul d'une distance tangentielle entre les images (Simard et al. (1998)). Le programme Matlab fait appel à un programme en C. L'intégration dans du code python plutôt que Matlab resterait à faire...
Les forêts aléatoires sont également une approche raisonnable, à moindre coût de développement, sur ces données. Analyser en détail la liste des paramètres proposés dans l'implémentation de cet algorithme. Consulter pour ce faire la documentation en ligne.
Les valeurs par défaut des paramètres sont utilisées sauf pour le nombre d'arbres: 100 au lieu de 10, et le nombre de processeurs utilisés: -1 au lieu de 1 (tous sont utilisés sauf 1 pour le système). Attention, tous les paramètres disponibles ne sont pas listés.
In [ ]:
from sklearn.ensemble import RandomForestClassifier
tps0 = time.clock()
rf = RandomForestClassifier(n_estimators=100,
criterion='gini', max_depth=None, min_samples_split=2,
min_samples_leaf=1, max_features='auto', max_leaf_nodes=None,
bootstrap=True, oob_score=True, n_jobs=-1,random_state=None, verbose=0)
rf.fit(Dtrain,Ltrain)
tps1 = time.clock()
print("Temps de configutration RF :" ,tps1 - tps0)
In [ ]:
# erreur out-of-bag
erreur_oob=1-rf.oob_score_
tps2 = time.clock()
print("Temps execution RF :", tps2 - tps0, "Erreur oob:", erreur_oob)
In [ ]:
# erreur sur l'échantillon test
1-rf.score(Dtest,Ltest)
In [ ]:
cm = confusion_matrix(Ltest, rf.predict(Dtest))
print(cm)
Comme pour les $k$ plus proches voisins, il serait utile d'optimiser certains paramètres dont le nombre d'arbres et sans doute max_features. L'optimisation de l'erreur out-of-bag plutôt qu'une procédure lourde de validaiton croisée serait bienvenue. D'autre part, la restriction de la profondeur max des arbres pourrait réduire sensiblement les temps de calcul mais cela ne semble pas nécessaire d'autant que c'est un paramètre critique pour la qualité de la prévision.
Le taux d'erreur de 3% obtenu sans effort d'optimisation est tout à fait correct au regard du temps passé en développement ! Plutôt que de chercher à l'optimiser, la suite du travail s'intéresse à l'effet de la taille de cet échantillon d'apprentissage sur la précision. La fonction learning_curve réalise ce calcul mais ne permet pas d'extraire le temps d'excution pour chaque taille. Une procédure rudimentaire est mise en oeuvre.
In [ ]:
from sklearn.cross_validation import train_test_split
# tailles croissantes de l'échantillon d'apprentissage
arrayErreur=np.empty((12,3))
nArbres=100
for i in range(1,13):
n=5000*i
arrayErreur[i-1,0]=n
if i==12:
n=59999
Xtrain,Xdrop,ytrain,ydrop=train_test_split(Dtrain,Ltrain,train_size=n)
tps1 = time.clock()
rf = RandomForestClassifier(n_estimators=nArbres,
criterion='gini', max_depth=None, min_samples_split=2,
min_samples_leaf=1, max_features='auto', max_leaf_nodes=None,
bootstrap=True, oob_score=True, n_jobs=-1,random_state=None, verbose=0)
rf.fit(Xtrain,ytrain)
tps2=time.clock()
arrayErreur[i-1,2]=1-rf.score(Dtest,Ltest)
arrayErreur[i-1,1]=tps2 - tps1
dataframeErreur1=pd.DataFrame(arrayErreur,columns=["Taille","Temps","Erreur"])
print(dataframeErreur1)
In [ ]:
# Graphes superposés
from __future__ import division
from scipy import *
from pylab import *
x = linspace(5,60,12)
fig = plt.figure()
# premier graphe
ax1 = fig.add_subplot(111)
ax1.plot(x,dataframeErreur1["Temps"] , '-b', label=ur"Temps",lw=1.5)
# absisses communes
xlim(0,65)
xlabel(ur"Taille échantillon x 1000", color='b', fontsize=16)
ylim(0, 70)
ylabel(ur"Secondes", color='b', fontsize=16)
legend(loc=2)
# 2ème graphe
ax2 = ax1.twinx()
ax2.plot(x,dataframeErreur1["Erreur"] ,'--g', label=ur"Erreur",lw=1.5)
ylim(0, 0.1)
ylabel(ur"Taux d'erreur", color='g', fontsize=16)
legend(loc=1)
show()
In [ ]:
from sklearn.ensemble import RandomForestClassifier
from sklearn.cross_validation import train_test_split
# tailles croissantes de l'échantillon d'apprentissage
arrayErreur=np.empty((12,3))
nArbres=250
for i in range(1,13):
n=5000*i
arrayErreur[i-1,0]=n
if i==12:
n=59999
Xtrain,Xdrop,ytrain,ydrop=train_test_split(Dtrain,Ltrain,train_size=n)
tps1 = time.clock()
rf = RandomForestClassifier(n_estimators=nArbres,
criterion='gini', max_depth=None, min_samples_split=2,
min_samples_leaf=1, max_features='auto', max_leaf_nodes=None,
bootstrap=True, oob_score=True, n_jobs=-1,random_state=None, verbose=0)
rf.fit(Xtrain,ytrain)
tps2=time.clock()
arrayErreur[i-1,2]=1-rf.score(Dtest,Ltest)
arrayErreur[i-1,1]=tps2 - tps1
dataframeErreur=pd.DataFrame(arrayErreur,columns=["Taille","Temps","Erreur"])
print(dataframeErreur)
In [ ]:
# Graphes supersosés
from __future__ import division
from scipy import *
from pylab import *
x = linspace(5,60,12)
fig = plt.figure()
# premier graphe
ax1 = fig.add_subplot(111)
ax1.plot(x,dataframeErreur["Temps"] , '-b', label=ur"Temps",lw=1.5,marker=".",markersize=6)
# absisses communes
xlim(0,65)
xlabel(ur"Taille échantillon x1000", fontsize=15)
ylim(0, 100)
ylabel(ur"Temps (s)", color='b', fontsize=15)
legend(loc=2)
# 2ème graphe
ax2 = ax1.twinx()
ax2.plot(x,dataframeErreur["Erreur"] ,'--',color='black', label=ur"Erreur",lw=1.5,marker=".",markersize=6)
ylim(0, 0.1)
ylabel(ur"Erreur (%)", fontsize=15)
legend(loc=1)
show()
Comparer les résultats obtenus (temps, précision) avec R.