Travaux pratiques - Comparaison de classifieurs dans Spark ML

Références externes utiles :

L'objectif de cette séance de TP est de comparer, sur un même jeu de données, quatre classifieurs disponibles dans Spark ML : la machine à vecteurs de support linéaire, la régression logistique, la forêt aléatoire et le perceptron multicouche (réseau de neurones). Pour chaque méthode nous construisons un pipeline, cherchons les valeurs optimales des hyperparamètres par validation croisée (avec grid search) et comparons les performances de généralisation sur un ensemble de test commun. La méthodologie est celle vue en cours ; il est utile de relire la section sur les pipelines et la sélection de modèle.

Nous travaillons sur le jeu de données Covertype (couverture forestière), issu de l'archive UCI. Il contient 581 012 observations décrites par 54 variables (10 quantitatives et 44 binaires issues du codage disjonctif de variables catégorielles), à classer en 7 types de couverture forestière.

Mise en place

Pour préparer les répertoires et télécharger les données, ouvrez un terminal et entrez :

%%bash
mkdir -p tpApprentissageLargeEchelle/data
wget -nc https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass/covtype.bz2 -P tpApprentissageLargeEchelle/data/
bzip2 -d tpApprentissageLargeEchelle/data/covtype.bz2
pyspark

Lecture et préparation des données

# Lecture des données en format libsvm
donnees = spark.read.format("libsvm") \
               .option("numFeatures", "54") \
               .load("tpApprentissageLargeEchelle/data/covtype")
donnees.printSchema()
donnees.show(3, truncate=False)
# Distribution des classes
donnees.groupBy("label").count().orderBy("label").show()

Question :

Combien d'observations contient chaque classe ? La distribution est-elle équilibrée ? Qu'est-ce que cela implique pour l'évaluation des modèles ?

Correction :

Le jeu est fortement déséquilibré, les classes 1 et 2 représentant ensemble env. 85 % des observations. Un taux de bonne classification global peut donc être trompeur car un classifieur qui prédit systématiquement la classe majoritaire obtiendrait un score élevé. L'AUC (aire sous la courbe ROC) ou la précision par classe sont des mesures plus adaptées dans ce contexte.

Nous réduisons le nombre d'observations par échantillonnage stratifié pour maintenir des temps de calcul raisonnable et équilibrer la distribution des classes.

from pyspark.sql import functions as F

# Taux d'échantillonnage par classe (environ 10 %)
#fractions = {float(i): 0.10 for i in range(1, 8)}    # réduction,sans équilibrage
#fractions = {1: 0.02, 2: 0.02, 3: 0.14, 4: 1, 5: 0.5, 6: 0.3, 7: 0.2}  # équilibrage, réduction
fractions = {1: 0.01, 2: 0.01, 3: 0.07, 4: 0.9, 5: 0.25, 6: 0.15, 7: 0.1}  # équlibrage, plus forte réduction

donnees_ech = donnees.stat.sampleBy("label", fractions, seed=42).cache()

print(f"Taille de l'échantillon : {donnees_ech.count()} observations")
donnees_ech.groupBy("label").count().orderBy("label").show()

Note :

Les étiquettes de classe dans Covertype sont entre 1 et 7, or Spark ML exige des étiquettes dans \(\{0, \ldots, K-1\}\), il faut donc les décaler de 1 :

donnees_ech = donnees_ech.withColumn("label", F.col("label") - 1)

Nous divisons ensuite le jeu en ensemble d'apprentissage (70 %) et ensemble de test (30 %). L'ensemble de test est mis de côté et ne sera utilisé qu'à la fin pour l'évaluation finale de chaque modèle.

partitions = donnees_ech.randomSplit([0.7, 0.3], seed=100)
apprentissage = partitions[0].cache()
test = partitions[1].cache()
print(f"Apprentissage : {apprentissage.count()} | Test : {test.count()}")

Standardisation des variables

Les méthodes linéaires (SVM, régression logistique) sont sensibles à l'échelle des variables. Nous définissons un StandardScaler qui sera inclus dans tous les pipelines :

from pyspark.ml.feature import StandardScaler

scaler = StandardScaler(inputCol="features",
                        outputCol="scaledFeatures",
                        withStd=True, withMean=True)

Question :

Pourquoi applique-t-on la standardisation à l'intérieur du pipeline plutôt qu'avant de séparer apprentissage et test ? Quel problème introduirait une standardisation calculée sur la totalité des données ?

Correction :

En calculant la standardisation sur la totalité des données (apprentissage + test), on introduit une fuite de données (data leakage) : les paramètres du scaler (moyennes et écarts-types) sont calculés en utilisant des informations issues des données de test, qui ne sont censées être « vues » qu'au moment de l'évaluation. Le modèle peut alors être optimisé implicitement pour ces données, ce qui fausse l'estimation des performances de généralisation.

Partie 1 : SVM linéaire

from pyspark.ml import Pipeline
from pyspark.ml.classification import LinearSVC, OneVsRest
from pyspark.ml.tuning import CrossValidator, ParamGridBuilder
from pyspark.ml.evaluation import MulticlassClassificationEvaluator
import numpy as np

Spark ML implémente LinearSVC pour la discrimination binaire uniquement. Pour le problème multiclasse à 7 classes, nous utilisons le schéma de décomposition un contre tous (One-vs-Rest, OvR) : 7 classifieurs binaires sont appris en parallèle, chacun distinguant une classe de toutes les autres. La classe retenue est celle dont le classifieur produit la valeur de décision la plus élevée.

svm = LinearSVC(featuresCol="scaledFeatures", labelCol="label",
                maxIter=20)
ovr_svm = OneVsRest(classifier=svm, labelCol="label",
                    featuresCol="scaledFeatures",
                    predictionCol="prediction")

pipeline_svm = Pipeline(stages=[scaler, ovr_svm])

# Grille : valeurs de la constante de régularisation C
# (regParam correspond à l'inverse de C pour LinearSVC dans Spark)
grid_svm = ParamGridBuilder() \
               .addGrid(svm.regParam, [0.01, 0.1]) \
               .build()

evaluateur = MulticlassClassificationEvaluator(
                 labelCol="label", predictionCol="prediction",
                 metricName="accuracy")

cv_svm = CrossValidator(estimator=pipeline_svm,
                        estimatorParamMaps=grid_svm,
                        evaluator=evaluateur,
                        numFolds=5,
                        seed=42)

Question :

Combien de classifieurs binaires sont appris au total pour cette étape de grid search avec validation croisée 5-fold ?

Correction :

La grille contient 3 valeurs pour regParam et la validation croisée 5-fold entraîne 5 modèles par point de grille. Il y a donc au total \(3 \times 5 = 15\) modèles. Pour OvR avec 7 classes, chaque modèle implique 7 classifieurs binaires, soit \(15 \times 7 = 105\) classifieurs binaires au total.

import time
t0 = time.time()
modele_svm = cv_svm.fit(apprentissage)
t_svm = time.time() - t0
print(f"SVM - temps d'entraînement : {t_svm:.1f} s")
print(f"Meilleur regParam : "
      f"{modele_svm.getEstimatorParamMaps()[np.argmax(modele_svm.avgMetrics)][svm.regParam]}")

acc_svm = evaluateur.evaluate(modele_svm.transform(test))
print(f"SVM accuracy sur test : {acc_svm:.4f}")

Partie 2 : Régression logistique

La régression logistique multinomiale (softmax) est directement disponible pour les problèmes multiclasses dans Spark ML, sans schéma OvR.

from pyspark.ml.classification import LogisticRegression

lr = LogisticRegression(featuresCol="scaledFeatures", labelCol="label",
                        maxIter=50, family="multinomial")

pipeline_lr = Pipeline(stages=[scaler, lr])

grid_lr = ParamGridBuilder() \
              .addGrid(lr.regParam,    [0.001, 0.01, 0.1]) \
              .addGrid(lr.elasticNetParam, [0.0, 0.5]) \
              .build()

cv_lr = CrossValidator(estimator=pipeline_lr,
                       estimatorParamMaps=grid_lr,
                       evaluator=evaluateur,
                       numFolds=5,
                       seed=42)

Question :

Le paramètre elasticNetParam de la régression logistique varie entre 0 et 1. Que représente-t-il ? A quoi correspondent les valeurs particulières 0,0 et 1,0 ?

Correction :

elasticNetParam (noté \(\alpha\)) contrôle le mélange des pénalisations L1 (LASSO) et L2 (ridge) dans le terme de régularisation : la pénalité est \(\alpha \|\mathbf{w}\|_1 + \frac{1-\alpha}{2}\|\mathbf{w}\|_2^2\). Pour \(\alpha = 0\) on obtient une régularisation L2 pure (ridge), les coefficients tendent vers 0 sans être exactement nuls. Pour \(\alpha = 1\) on obtient une régularisation L1 pure (LASSO) qui favorise la parcimonie : certains coefficients deviennent exactement nuls. Des valeurs intermédiaires combinent les deux effets (elastic net).

t0 = time.time()
modele_lr = cv_lr.fit(apprentissage)
t_lr = time.time() - t0
print(f"Régression logistique - temps d'entraînement : {t_lr:.1f} s")

best_lr = modele_lr.getEstimatorParamMaps()[np.argmax(modele_lr.avgMetrics)]
print(f"Meilleur regParam : {best_lr[lr.regParam]}, "
      f"elasticNetParam : {best_lr[lr.elasticNetParam]}")

acc_lr = evaluateur.evaluate(modele_lr.transform(test))
print(f"Régression logistique accuracy sur test : {acc_lr:.4f}")

Partie 3 : Forêt aléatoire

La forêt aléatoire est une méthode d'ensemble qui agrège les prédictions d'un grand nombre d'arbres de décision appris sur des sous-échantillons des données et des sous-ensembles de variables. Les arbres de décision ne nécessitent pas de standardisation préalable.

from pyspark.ml.classification import RandomForestClassifier

rf = RandomForestClassifier(featuresCol="features", labelCol="label",
                            numTrees=50, seed=42)

# Pas de standardisation nécessaire pour la forêt aléatoire
pipeline_rf = Pipeline(stages=[rf])

grid_rf = ParamGridBuilder() \
              .addGrid(rf.maxDepth,  [5, 10]) \
              .addGrid(rf.numTrees,  [50, 100]) \
              .build()

cv_rf = CrossValidator(estimator=pipeline_rf,
                       estimatorParamMaps=grid_rf,
                       evaluator=evaluateur,
                       numFolds=5,
                       seed=42)

Question :

Quel est l'impact des hyperparamètres maxDepth et numTrees sur la capacité de généralisation du modèle ?

Correction :

maxDepth contrôle la profondeur maximale de chaque arbre. Un arbre peu profond peut être biaisé (sous-apprentissage), un arbre très profond peut sur-apprendre sur les données d'apprentissage. numTrees contrôle le nombre d'arbres dans la forêt : plus il y en a, plus la variance de la prédiction est réduite, mais au prix d'un temps de calcul plus élevé. En pratique, sa valeur se stabilise assez vite avec le nombre d'arbres.

t0 = time.time()
modele_rf = cv_rf.fit(apprentissage)
t_rf = time.time() - t0
print(f"Forêt aléatoire - temps d'entraînement : {t_rf:.1f} s")

best_rf = modele_rf.getEstimatorParamMaps()[np.argmax(modele_rf.avgMetrics)]
print(f"Meilleur maxDepth : {best_rf[rf.maxDepth]}, "
      f"numTrees : {best_rf[rf.numTrees]}")

acc_rf = evaluateur.evaluate(modele_rf.transform(test))
print(f"Forêt aléatoire accuracy sur test : {acc_rf:.4f}")

Partie 4 : Réseau de neurones (MLP)

Spark ML propose un perceptron multicouche (Multilayer Perceptron Classifier, MLP). L'architecture du réseau est définie par le vecteur layers : le premier élément est la dimension d'entrée (54 variables), le dernier est le nombre de classes (7), les éléments intermédiaires sont les dimensions des couches cachées. Il n'y a pas d'hyperparamètre explicite de régularisation (par ex. weight decay) et si on souhaite employer early stopping il faut le coder soi-même.

Attention, le temps de calcul est beaucoup plus long ici, il est envisageable d'éviter cet apprentissage durant le TP.

from pyspark.ml.classification import MultilayerPerceptronClassifier

mlp = MultilayerPerceptronClassifier(
          featuresCol="scaledFeatures", labelCol="label",
          layers=[54, 64, 32, 7],
          maxIter=100, seed=42)

pipeline_mlp = Pipeline(stages=[scaler, mlp])

grid_mlp = ParamGridBuilder() \
               .addGrid(mlp.layers, [[54, 64, 32, 7],
                                     [54, 128, 64, 7]]) \
               .addGrid(mlp.stepSize, [0.01, 0.05]) \
               .build()

cv_mlp = CrossValidator(estimator=pipeline_mlp,
                        estimatorParamMaps=grid_mlp,
                        evaluator=evaluateur,
                        numFolds=5,
                        seed=42)

Question :

Calculez le nombre de paramètres (poids et biais) de l'architecture [54, 64, 32, 7]. Comparez à celui de la régression logistique multinomiale pour le même problème. Qu'en déduisez-vous sur le risque de sur-apprentissage et sur le volume de données nécessaire ?

Correction :

Pour [54, 64, 32, 7] : la couche 1 a \(54 \times 64 + 64 = 3520\) paramètres ; la couche 2 a \(64 \times 32 + 32 = 2080\) ; la couche 3 a \(32 \times 7 + 7 = 231\) ; au total il y donc 5831 paramètres. Pour la régression logistique multinomiale nous avions \(54 \times 7 + 7 = 385\) paramètres. Le MLP a environ 15 fois plus de paramètres à estimer, ce qui augmente le risque de sur-apprentissage et exige davantage de données d'apprentissage (ou une régularisation plus forte) pour obtenir une bonne généralisation.

t0 = time.time()
modele_mlp = cv_mlp.fit(apprentissage)
t_mlp = time.time() - t0
print(f"MLP - temps d'entraînement : {t_mlp:.1f} s")

best_mlp = modele_mlp.getEstimatorParamMaps()[np.argmax(modele_mlp.avgMetrics)]
print(f"Meilleur layers : {best_mlp[mlp.layers]}, "
      f"stepSize : {best_mlp[mlp.stepSize]}")

acc_mlp = evaluateur.evaluate(modele_mlp.transform(test))
print(f"MLP accuracy sur test : {acc_mlp:.4f}")

Partie 5 : Comparaison et analyse

Tableau récapitulatif

resultats = {
    "SVM (OvR)":             (acc_svm,  t_svm),
    "Régression logistique": (acc_lr,   t_lr),
    "Forêt aléatoire":       (acc_rf,   t_rf),
    "MLP":                   (acc_mlp,  t_mlp),
}

print(f"{'Méthode':<25} {'Accuracy (test)':>16} {'Temps entr. (s)':>16}")
print("-" * 58)
for nom, (acc, t) in sorted(resultats.items(), key=lambda x: -x[1][0]):
    print(f"{nom:<25} {acc:>16.4f} {t:>16.1f}")

Analyse par classe

Les classes présentes sont :
  • Classe 0 : Spruce/Fir (épicea/sapin)

  • Classe 1 : Lodgepole Pine (pin tordu)

  • Classe 2 : Ponderosa Pine (pin ponderosa)

  • Classe 3 : Cottonwood/Willow (peuplier/saule))

  • Classe 4 : Aspen (peuplier faux-tremble))

  • Classe 5 : Douglas-fir (Sapin Douglas)

  • Classe 6 : Krummholz (arbres nains)

Un taux de bonne classification global masque des différences importantes selon les classes, surtout sur un jeu déséquilibré. Nous calculons ici la précision et le rappel par classe pour la meilleure méthode.

# À remplacer par le meilleur modèle selon les résultats obtenus
meilleur_modele = modele_rf    # à adapter

predictions = meilleur_modele.transform(test)

# Évaluation par classe
for metrique in ["precisionByLabel", "recallByLabel", "fMeasureByLabel"]:
    eval_mc = MulticlassClassificationEvaluator(
                  labelCol="label", predictionCol="prediction",
                  metricName=metrique)
    valeurs = eval_mc.evaluate(predictions)
    print(f"{metrique} : {valeurs:.4f}")
# Matrice de confusion
from pyspark.mllib.evaluation import MulticlassMetrics

preds_et_labels = predictions.select("prediction", "label") \
                              .rdd.map(lambda r: (r.prediction, float(r.label)))
metrics = MulticlassMetrics(preds_et_labels)
print("Matrice de confusion :")
print(metrics.confusionMatrix().toArray().astype(int))

Question :

Quelles sont les classes les plus difficiles à prédire ? Ces difficultés sont-elles symétriques (deux classes se confondent mutuellement) ou asymétriques ? Quelle hypothèse cela suggère-t-il sur la structure des données pour ces classes ?

Correction :

Sans équilibrage, les classes les moins représentées (4, 5, 6 et 7) sont généralement les plus difficiles à prédire. Les confusions sont plus fréquentes entre des types de couverture forestière proches (par ex. épicéa vs sapin dans des conditions géographiques similaires). Des confusions asymétriques entre deux classes A et B (A souvent prédit comme B mais pas l'inverse) suggèrent que A est un sous-ensemble ou une variante moins typique de B, ce qui est cohérent avec le fait que des types forestiers voisins partagent des caractéristiques géographiques similaires.

Question (bilan) :

Comparez les quatre méthodes selon trois critères : accuracy sur les données de test, temps d'entraînement et interprétabilité du modèle. Y a-t-il une méthode qui domine sur tous les critères ? Comment le choix dépend-il des contraintes applicatives (contrainte de temps réel, besoin d'explication des décisions, volume de données disponible) ?

Correction :

En général sur Covertype la forêt aléatoire et le MLP obtiennent les meilleures accuracies, devant la régression logistique et la SVM linéaire. La forêt aléatoire offre en plus une mesure d'importance des variables utile pour l'interprétation. Le MLP est souvent le plus lent à entraîner. La SVM et la régression logistique sont les plus rapides lors de l'inférence (complexité linéaire en \(p\)).

Aucune méthode ne domine sur tous les critères : la forêt aléatoire offre un bon compromis accuracy/interprétabilité (importance des variables), la régression logistique est la plus simple à mettre en production, le MLP peut capturer des non-linéarités (comme la forêt aléatoire, d'ailleurs) mais est une boîte noire. Le choix dépend du contexte : si contrainte légale d'explication alors régression logistique ou forêt aléatoire ; si maximisation de la performance alors MLP ou forêt aléatoire ; si contrainte de faible latence à l'inférence alors SVM ou régression logistique.

Importance des variables (forêt aléatoire)

La forêt aléatoire fournit une mesure d'importance des variables qui est utile pour comprendre quelles variables contribuent le plus à la discrimination.

# Récupération du modèle RandomForest depuis le meilleur pipeline
rf_model = modele_rf.bestModel.stages[-1]
importances = rf_model.featureImportances.toArray()

# Affichage des 10 variables les plus importantes
top10 = np.argsort(importances)[::-1][:10]
print("Top 10 variables les plus importantes :")
for rang, idx in enumerate(top10, 1):
    print(f"  {rang:2d}. variable {idx:3d} : importance = {importances[idx]:.4f}")

Question :

Les variables les plus importantes sont-elles parmi les 10 variables quantitatives initiales (indices 0 à 9) ou parmi les variables binaires issues du codage disjonctif (indices 10 à 53) ? Qu'est-ce que cela suggère sur l'utilité respective des deux types de variables pour ce problème ?

Correction :

Dans Covertype les variables quantitatives (élévation, pente, distance aux points d'eau, ensoleillement) ont en général une importance plus élevée que les variables binaires issues du codage disjonctif des types de sol. L'élévation est typiquement la variable la plus discriminante car les types de couverture forestière sont fortement corrélés à l'altitude. Cela indique que les variables continues originales portent l'essentiel de l'information discriminanteet les variables binaires une information complémentaire.