Travaux pratiques - Evaluation et sélection de modèles décisionnels

L’objectif de cette séance de travaux pratiques est de montrer l’utilisation des techniques de validation croisée pour l’évaluation et la comparaison de modèles décisionnels, ainsi que des méthodes de recherche de valeurs pour les hyper-paramètres (comme le coefficient de régularisation).

Références externes utiles :

Estimation des performances par validation croisée

Afin d’illustrer l’utilisation de la validation croisée, nous considérons un problème de classement similaire à celui examiné lors de la précédente séance. Nous considérons un jeu de données pour lequel nous disposons de l’information de supervision (les étiquettes des classes) et nous le partitionnons en un ensemble d’apprentissage et un ensemble de test. Nous utiliserons les Iris de Fisher comme jeu de données. Pour rappel, celui-ci consiste à identifier à quelle espèce appartient une plante parmi trois possibilités (Iris Setosa, Iris Virginica et Iris Versicolor) à partir de la longueur et la largeur des sépales et des pétales.

Nous employons des perceptrons multi-couches (PMC) avec une seule couche cachée de 100 neurones et une valeur \(\alpha = 1\) pour la constante de régularisation (pondération du terme d’oubli ou weight decay). Il n’est pas indispensable pour cette séance de savoir en détails comment fonctionne ces réseaux de neurones, ces modèles décisionnels seront revus plus tard dans le cours.

La validation croisée sera utilisée pour estimer les performances de généralisation à partir de l’ensemble d’apprentissage. Ensuite, cette estimation sera comparée à l’estimation obtenue sur l’ensemble de test que nous avions mis de côté au départ. Les explications sur la validation croisée et sa mise en œuvre dans Scikit-learn se trouvent dans la documentation.

# importations
import numpy as np
import matplotlib.pyplot as plt
import sklearn

# chargement du jeu de données Iris
from sklearn import datasets
iris = datasets.load_iris()
data, labels = iris.data, iris.target

# découpage initial en données d'apprentissage et données de test
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(data, labels, test_size=0.5)

La fonction train_test_split de scikit-learn nous permet de diviser aléatoirement le jeu de données en deux partitions train (apprentissage) et test (évaluation) selon des proportions arbitraires.

Comme d’habitude, il est intéressant de visualiser les données à notre disposition. Nous pouvons construire les nuages de points en deux dimensions à l’aide de Matplotlib :

fig = plt.figure(figsize=(16, 12))
n_features = data.shape[-1]
n_plots = 6
idx = 1
cmp = np.array(['r', 'g', 'b'])
for dim1 in range(0, n_features):
    for dim2 in range(dim1+1, n_features):
        fig.add_subplot(2, n_plots // 2, idx)
        plt.scatter(X_train[:, dim1], X_train[:, dim2],c=cmp[y_train], s=50, edgecolors='none')
        plt.scatter(X_test[:,  dim1], X_test[:, dim2], c='none',  s=50, edgecolors=cmp[y_test])
        plt.xlabel(iris.feature_names[dim1])
        plt.ylabel(iris.feature_names[dim2])
        idx += 1
plt.show()

Question

Pourquoi a-t-on six nuages de points différents ? Que représente chaque figure ?

Correction

Il y a un nuage de points pour chaque paire de variables. Chaque figure représente donc la projection du jeu de données selon deux dimensions.

Pour ce problème de classement, nous allons utiliser un PMC avec les valeurs par défaut proposés par scikit-learn pour le nombre de couches (1) et de neurones (100).

# emploi de PMC
from sklearn.neural_network import MLPClassifier
clf = MLPClassifier(solver='lbfgs', alpha=1, tol=5e-3)

Pour estimer l’erreur de généralisation, nous allons utiliser la validation croisée K-fold. Cela va nous permettre d’ajuster, si besoin, les hyperparamètres du modèle décisionnel. scikit-learn implémente diverses stratégies de validation croisée dans le module sklearn.model_selection. Commençons par expérimenter avec l’approche K-fold. L’objet KFold dispose d’une méthode .split() qui permet de générer les listes des indices des observations à utiliser pour le sous-jeu d’apprentissage et pour le sous-jeu de validation. Plus de détails sur cet objet se trouvent dans la documentation de K-Fold.

# KFold pour différentes valeurs de k
from sklearn.model_selection import KFold

# valeurs de k
n_folds = np.array([2, 3, 5, 7, 10, 13, 16, 20])

# préparation des listes pour stocker les résultats
cv_scores = []
cv_scores_std = []

for k in n_folds:    # pour chaque valeur de k
    kf = KFold(n_splits=k)
    scores = []
    # apprentissage puis évaluation d'un modèle sur chaque split
    for train_idx, val_idx in kf.split(X_train):
        # apprentissage avec .fit()
        clf.fit(X_train[train_idx], y_train[train_idx])
        scores.append(clf.score(X_train[val_idx], y_train[val_idx]))
    # calcul de la moyenne et de l'écart-type des performances obtenues
    cv_scores.append(np.mean(scores))
    cv_scores_std.append(np.std(scores))

cv_scores, cv_scores_std = np.array(cv_scores), np.array(cv_scores_std)

# affichage performance moyenne +- 1 écart-type pour chaque k
plt.figure(figsize=(8, 6))
plt.plot(n_folds, cv_scores, 'b')
plt.fill_between(n_folds, cv_scores + cv_scores_std, cv_scores - cv_scores_std, alpha=0.5)
plt.xlabel("Valeur de $k$ du K-Fold")
plt.ylabel("Score moyen")
plt.xlim(2, max(n_folds))
plt.xticks(n_folds)
plt.title("Erreur de généralisation estimée en fonction de $k$")
plt.show()

Question :

Que constatez-vous en examinant ce graphique ? Ajoutez des valeurs pour k (par ex. 40, 100, attention ce sera plus long…) et examinez de nouveau le graphique.

Correction :

On constate que, lors de l’augmentation de k, la performance moyenne se stabilise mais la variance augmente. Cela s’explique par le fait que, lorsque la valeur de k augmente, l’évaluation est faite (c’est à dire la moyenne de l’erreur est calculée) sur de moins en moins de données. La variance augmente encore pour des valeurs supérieures de k.

Question :

Pour chaque modèle appris par validation croisée k-fold, ajoutez son évaluation sur les données de test mises de côté au départ X_test, y_test. Affichez les courbes sur le même graphique. Que constatez-vous ?

Correction :

On ajoute des listes pour stocker ces résultats et on affiche leurs contenus :

kcvscores = list()
kcvscores_std = list()
testscores = list()
testscores_std = list()

kcvfs = [2, 3, 4, 5, 10, 30]

for kcvf in kcvfs:    # pour chaque valeur de k
  kf = KFold(n_splits=kcvf)
  these_scores = list()
  these_test_scores = list()
  # apprentissage puis évaluation d'un modèle sur chaque split
  for train_idx, test_idx in kf.split(X_train):
    clf.fit(X_train[train_idx], y_train[train_idx])
    these_scores.append(clf.score(X_train[test_idx], y_train[test_idx]))
    these_test_scores.append(clf.score(X_test, y_test))
  # calcul de la moyenne et de l'écart-type des performances obtenues
  kcvscores.append(np.mean(these_scores))
  kcvscores_std.append(np.std(these_scores))
  testscores.append(np.mean(these_test_scores))
  testscores_std.append(np.std(these_test_scores))

# création de np.array à partir des listes
kcvscores, kcvscores_std = np.array(kcvscores), np.array(kcvscores_std)
testscores, testscores_std = np.array(testscores), np.array(testscores_std)

# affichage performance moyenne +- 1 écart-type pour chaque k
plt.figure()
plt.plot(kcvfs, kcvscores, 'b')
plt.plot(kcvfs, kcvscores+kcvscores_std, 'b--')
plt.plot(kcvfs, kcvscores-kcvscores_std, 'b--')
plt.plot(kcvfs, testscores, 'g')
plt.plot(kcvfs, testscores+testscores_std, 'g--')
plt.plot(kcvfs, testscores-testscores_std, 'g--')

Les résultats montrent que l’estimation de l’erreur de généralisation par validation croisée sur les données d’apprentissage (courbes en bleu) reste en général optimiste par rapport à l’estimation sur des données de test supplémentaires (courbes en vert). Aussi, la variance des estimations sur les données de test est comparativement faible car ces données sont ici aussi volumineuses que les données d’apprentissage (test_size=0.5).

Question :

Réalisez l’estimation des performances en utilisant la validation croisée leave one out (LOO). Que constatez-vous en comparant les résultats de k-fold et de leave one out ?

Correction :

# LOO
from sklearn.model_selection import LeaveOneOut
loo = LeaveOneOut()
loo.get_n_splits(X_train)
these_scores = list()
for train_idx, test_idx in loo.split(X_train):
  clf.fit(X_train[train_idx], y_train[train_idx])
  these_scores.append(clf.score(X_train[test_idx], y_train[test_idx]))

np.mean(these_scores)
np.std(these_scores)

On constate que l’écart-type est bien plus élevé pour l’estimation leave one out que pour les estimations k-fold (pour toutes les valeurs considérées ici pour k).

Recherche des meilleures valeurs pour les hyperparamètres

Une des problématiques centrales en modélisation décisionnelle est de pouvoir déterminer les hyperparamètres du modèle qui débouchent sur les meilleurs performances. Formellement, cela revient à chercher la famille paramétrique qui obtient l’erreur de généralisation la plus faible. Néanmoins, comme nous l’avons vu, nous ne pouvons pas utiliser le jeu d’évaluation pour comparer des centaines de modèles : l’estimation de l’erreur généralisation deviendrait bien trop optimiste. Nous allons donc utiliser la validation croisée pour comparer nos modèles, puis seulement une fois les meilleures valeurs des hyperparamètres déterminées pourrons nous estimer l’erreur de généralisation sur le jeu de test.

La méthode la plus naïve pour la recherche d’hyperparamètres est la recherche systématique, ou grid search. Elle consiste à explorer toutes les combinaisons de valeurs possibles dans un ensemble donné.

Utilisons grid search pour trouver les meilleures valeurs de deux hyperparamètres pour les PMC dans la même tâche de classement que précédemment. Ces hyperparamètres sont le nombre de neurones dans l’unique couche cachée du PMC et la valeur de la constante de régularisation (par weight decay), \(\alpha\).

# importations
import numpy as np
import matplotlib.pyplot as plt

# chargement des données iris
from sklearn import datasets
data, labels = datasets.load_iris(return_X_y=True)

# découpage initial en données d'apprentissage et données de test
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(data, labels, test_size=0.5)

# emploi de PMC
from sklearn.neural_network import MLPClassifier

Afin d’utiliser la recherche dans une grille et la validation croisée pour comparer les modèles obtenus avec toutes les combinaisons de valeurs pour les hyperparamètres, il est nécessaire d’employer GridSearchCV. Cet estimateur scikit-learn combine automatiquement une rechercher par grille avec une validation croisée K-fold.

from sklearn.model_selection import GridSearchCV

Il est nécessaire d’indiquer dans un « dictionnaire » quels sont les hyperparamètres dont on souhaite explorer les valeurs et quelles sont les différentes valeurs à évaluer. Chaque entrée du dictionnaire consiste en une chaîne de caractères qui contient le nom de l’hyperparamètre tel qu’il est défini dans l’estimateur employé. Nous nous servirons ici de MLPClassifier, les noms des paramètres peuvent donc être trouvés dans la présentation de cette classe. Nous considérons ici seulement deux paramètres, hidden_layer_sizes (nombre de neurones dans l’unique couche cachée) et alpha (la constante \(\alpha\) de régularisation par weight decay).

tuned_parameters = {'hidden_layer_sizes':[(5,), (20,), (50,), (100,), (150,), (200,)],
                    'alpha':   [0.001, 0.01, 1, 2]}

Dans l’appel de GridSearchCV nous indiquons ensuite pour MLPClassifier le solveur à utiliser systématiquement (qui n’est pas celui par défaut), ensuite le dictionnaire avec les valeurs des (hyper)paramètres à explorer et enfin le fait que c’est la validation croisée k-fold avec \(k=5\) qui est employée pour comparer les différents modèles.

clf = GridSearchCV(MLPClassifier(solver='lbfgs', tol=5e-3), tuned_parameters, cv=5)

# exécution de grid search
clf.fit(X_train, y_train)

Scikit-learn exécute alors le programme suivant :

  • à partir des listes de valeurs pour les différents (hyper)paramètres sont générées toutes les combinaisons de valeurs,

  • pour chaque combinaison, les performances des modèles correspondants sont évaluées par validation croisée 5-fold (appliquée uniquement sur les données d’apprentissage X_train, y_train),

  • sont sélectionnées les valeurs des (hyper)paramètres correspondant aux meilleures performances de validation croisée,

  • avec ces valeurs pour les (hyper)paramètres un nouvel apprentissage est réalisé avec la totalité des données de X_train, y_train (et non seulement \(\frac{k-1}{k}\) folds).

Les lignes suivantes permettent d’afficher les résultats : les paramètres du meilleur modèle avec clf.best_params_, ainsi que les résultats de validation croisée obtenus pour toutes les combinaisons de valeurs pour les (hyper)paramètres (clf.cv_results_ donne accès à ces informations et à bien d’autres).

print(clf.best_params_)

Nous pouvons également faire un affichage plus complet, traçant la surface du score moyen des différentes modèles en fonction de la combinaison des deux hyperparamètres :

n_hidden = np.array([s[0] for s in tuned_parameters['hidden_layer_sizes']])
alphas = np.array(tuned_parameters['alpha'])

# création de la grille des hyperparamètres
xx, yy = np.meshgrid(n_hidden, alphas)
scores = clf.cv_results_['mean_test_score'].reshape(xx.shape)

# affichage sous forme de fil de fer de la surface des résultats des modèles évalués
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8, 6))
ax = plt.axes(projection='3d')
ax.set_xlabel("Neurones cachés")
ax.set_ylabel("Régularisation $\\alpha$")
ax.set_zlabel("Taux de bon classement")
ax.plot_wireframe(xx, yy, scores)
plt.show()

Nous avons employé ici plot_wireframe car la lisibilité est meilleure qu’avec plot_surface.

Question :

Combien de PMC sont appris au total dans cet exemple ?

Correction :

Le nombre de combinaisons de (hyper)paramètres explorées est len(tuned_parameters['hidden_layer_sizes']) * len(tuned_parameters['alpha']) = 24. Pour chacune de ces combinaisons, cv=5 indique que \(k = 5\) PMC différents sont appris. Donc un total de 24 * 5 = 120.

Question :

Quelle est la signification du paramètre refit de GridSearchCV ?

Correction :

Comme indiqué dans la documentation, si la valeur de ce paramètre est True (valeur par défaut) alors, une fois trouvées les meilleures valeurs pour les hyperparamètres, un nouveau modèle est appris avec ces valeurs-là sur la totalité des \(N\) données d’apprentissage X_train, y_train (sans en exclure \(N/k\)). Ce modèle est directement accessible dans l’attribut .best_estimator_ et l’appel à .predict() sur l’instance de GridSearchCV (ici clf) permet de l’utiliser.

Question :

Examinez de façon plus complète le contenu de clf.cv_results_.

Correction :

Regarder les explications concernant cet attribut dans la documentation.

Question :

Evaluez le modèle sélectionné sur les données de test (X_test, y_test).

Correction :

Le paramètres refit étant par défaut True, le modèle appris avec les meilleures valeurs pour les hyperparamètres est directement accessible via l’instance de GridSearchCV (ici clf), donc pour l’évaluer sur les données de test il suffit d’écrire

clf.score(X_test, y_test)

Question :

L’aspect des résultats vous incite à affiner la grille ? Modifiez la grille et examinez les nouveaux résultats.

Correction :

Il est surtout intéressant d’affiner la grille autour des valeurs optimales pour les (hyper)paramètres, lues sur le graphique affiché ou obtenues avec clf.best_params_. Il faut définir une nouvelle grille plus fine autour de ce point et appeler de nouveau GridSearchCV.

Question :

Utilisez la recherche aléatoire avec RandomizedSerchCV. Le « budget » (nombre total de combinaisons évaluées) peut être fixé avec n_iter. Motivez le choix des lois employées pour le tirage des valeurs des deux (hyper)paramètres hidden_layer_sizes et alpha.

Correction :

Les distributions peuvent être choisies dans cette liste de scipy.stats. Les distributions continues doivent être préférées pour les paramètres continus (comme \(\alpha\) ici) et les distributions discrètes pour les paramètres discrets (comme le nombre de neurones dans la couche cachée). Les distributions uniformes (uniform, respectivement randint) sont la solution de facilité. Si des connaissances a priori nous permettent de préférer certains points de l’espace des paramètres, alors nous pouvons choisir d’autres distributions qui privilégient les voisinages de ces points. L’appel à RandomizedSerchCV aura la forme

rlf = RandomizedSerchCV(MLPClassifier(solver='lbfgs'), param_distributions=param_distrib, n_iter=50, cv=5)

param_distrib est le dictionnaire qui précise les distributions employées pour les différents (hyper)paramètres et n_iter=50 indique un « budget » de 50 essais.

Pour aller plus loin

D’autres stratégies plus avancées de sélection d’hyperparamètres sont envisageables. En effet, la recherche par grille, comme la recherche aléatoire, sont des méthodes relativement rudimentaires et plutôt coûteuses en ressources. Cela peut notamment poser problème lorsque les modèles décisionnels que l’on considère appartiennent à une famille paramétrique coûteuse en calculs, comme les réseaux de neurones profonds.

scikit-learn propose ainsi une implémentation expérimentale d’une recherche par dichotomie, permettant de concentrer les ressources sur les combinaisons d’hyperparamètres les plus prometteuses.

Question

Expérimenter avec HalvingGridSearchCV (cf. la documentation). Comparer notamment le score final et le temps de calcul à celui obtenu avec GridSearchCV.

Hors de scikit-learn, il existe également d’autres bibliothèques implémentant des stratégies avancées de sélection utilisant diverses heuristiques pour choisir les combinaisons d’hyperparamètres à explorer. Certaines de ces bibliothèques sont compatibles avec l’interface de scikit-learn, par exemple tune-sklearn ou encore sklearn-genetic-opt.