TP++ - Analyse factorielle discriminante - Étape descriptive

(les TP++ sont des compléments à l’enseignement de base et ne font pas partie des séances de travaux pratiques avec encadrement)

Références externes utiles :

L’objectif de cette séance de TP est de montrer une mise en œuvre avec Spark de l’étape descriptive de l’analyse factorielle discriminante (AFD) ou Linear Discriminant Analysis (LDA). Cette méthode d’analyse n’est pas (encore) proposée par Spark dans MLLib.

L’analyse factorielle discriminante traite des observations caractérisées par des variables quantitatives plus une variable nominale, de « classe », et cherche à trouver dans l’espace initial constitué par les variables quantitatives un sous-espace linéaire dans lequel la séparation entre classes est optimale. Une présentation rapide de l’AFD est dans cette section de RCP216 et une présentation plus détaillée dans cette section de RCP208.

L’approche employée ici pour réaliser l’AFD consiste à

  • effectuer une ACP normée des observations,

  • projeter les observations et les centres des classes sur les axes principaux correspondant à des valeurs propres non nulles,

  • appliquer aux observations une transformation qui égalise les variances des projections sur les axes principaux (whitening),

  • calculer les centres de gravité des classes,

  • appliquer une nouvelle ACP aux seuls centres de gravité, les axes principaux résultants étant ainsi les axes discriminants de l’analyse discriminante des observations,

  • projeter les observations sur les trois premiers axes discriminants.

La mise en œuvre de l’AFD est illustrée ici sur les données « textures », voir par exemple cette présentation. Ces données consistent en 5500 observations caractérisées par 40 variables numériques de départ et appartenant à 11 classes.

L’approche proposée ici peut être adaptée à la visualisation d’autres données vectorielles qui sont caractérisées aussi par des étiquettes de classes (adaptées aux modèles décisionnels de classement) ou des étiquettes de groupes (résultant de l’application de modèles descriptifs issus de la classification automatique, par ex. k-means).

La visualisation nécessite d’avoir au préalable importé la bibliothèque vegas-viz.

Cela peut se faire en ligne de commande lors du lancement de spark-shell : spark-shell --packages "org.vegas-viz:vegas_2.11:0.3.11","org.vegas-viz:vegas-spark_2.11:0.3.11".

Pour la lecture des données, il faut entrer les commandes suivantes dans une fenêtre terminal :

%%bash
mkdir -p tpafd/data
cd tpafd/data
wget http://cedric.cnam.fr/vertigo/Cours/RCP216/docs/texture.csv

Ensuite, si vous utilisez Zeppelin, il suffit d’exécuter la cellule ci-dessous. Sinon, dans une autre fenêtre terminal, vous devez vous positionner dans le répertoire ~/tpafd, lancer spark-shell et entrer dans spark-shell les commandes qui suivent (noter que la séquence d’opérations indiquée ci-dessous peut être améliorée) :

// Stage 1:
// - load the textures dataset
// - then find the centers of each of the 11 classes

import org.apache.spark.sql.types._

def genTexturesFields(from: Int, to: Int) = for (i <- from until to)
                      yield StructField("val"+i.toString, DoubleType, true)

val texturesFields = genTexturesFields(0, 40)

val texturesSchema = StructType(texturesFields).add("label", DoubleType, true)

// read the textures dataset in a DataFrame
val texturesDF = spark.read.format("csv").option("sep"," ").schema(texturesSchema)
                      .load("tpafd/data/texture.csv").cache()

// check reading (there are 500 samples for each of the 11 classes)
texturesDF.groupBy("label").count().show()

// define the aggregation operation by a Map
val colMap = {for (i <- 0 until 40) yield ("val"+i,"avg")}.toMap

// find the center of each of the 11 classes
val cgclassesDF = texturesDF.groupBy("label").agg(colMap)

// verification
cgclassesDF.groupBy("label").count().show

// Stage 2:
// - assemble the 1D features to prepare scaling and PCA,
// - then fit scaler, transform data, fit PCA
// - then project both data and class centers on the PCs

import org.apache.spark.ml.feature.VectorAssembler

// define assembler for data
val colsEntree = {for (i <- 0 until 40) yield "val"+i.toString}.toArray

val assembleur = new VectorAssembler().setInputCols(colsEntree).setOutputCol("features")

// assemble data
val texturesDFA = assembleur.transform(texturesDF).select("label","features")

// verification
texturesDFA.groupBy("label").count().show

// define assembler for class centers (1D features have different names)
val colsEntreeCG = {for (i <- 0 until 40) yield "avg(val"+i+")".toString}.toArray

val assembleurCG = new VectorAssembler().setInputCols(colsEntreeCG).setOutputCol("features")

//assemble class centers
val cgclassesDFA = assembleurCG.transform(cgclassesDF).select("label","features")

import org.apache.spark.ml.feature.StandardScaler
import org.apache.spark.ml.feature.PCA

// define scaler
val scaler = new StandardScaler().setInputCol("features")
                                 .setOutputCol("scaledFeatures")
                                 .setWithStd(true)
                                 .setWithMean(true)
// define PCA
val pca = new PCA().setInputCol("scaledFeatures").setOutputCol("pcaFeatures").setK(40)

// fit scaler and inspect
val scalerModel = scaler.fit(texturesDFA)

scalerModel.mean
scalerModel.std

// apply scaling to the data
val scaledTexturesDFA = scalerModel.transform(texturesDFA).select("label","scaledFeatures")
scaledTexturesDFA.printSchema()

// fit PCA and inspect
val pcaModel = pca.fit(scaledTexturesDFA)

pcaModel.explainedVariance

// Stage 3:
// - whiten the class centers,
// - then fit a PCA to the centers
//   (the resulting PCs are the discriminant components of the data)

import org.apache.spark.ml.feature.ElementwiseProduct
import org.apache.spark.ml.linalg.Vectors
import scala.math.sqrt

// prepare vector for PCA whitening
val wv = Vectors.dense(pcaModel.explainedVariance.toArray.map(x => 1/(sqrt(x))))

// preparer whitening transformer
val whitener = new ElementwiseProduct().setScalingVec(wv)
                                       .setInputCol("pcaFeatures")
                                       .setOutputCol("whitenedFeatures")

// apply scaling, PCA projection and then whitening to class centers
val cgclassesDFAwhitened = whitener.transform(pcaModel.transform(scalerModel
                                                      .transform(cgclassesDFA)))

// define new PCA on whitened class centers
//  the resulting PCs are the discriminant components of the initial data
val pcaCG = new PCA().setInputCol("whitenedFeatures").setOutputCol("ldaFeatures").setK(3)

// fit new PCA and inspect
val ldaModel = pcaCG.fit(cgclassesDFAwhitened)
ldaModel.explainedVariance

// Stage 4:
// - apply PCA1 projection, whitening and then PCA2 projection to the initial data
// - save resulting data to prepare visualization

// apply PCA1 projection, whitening and then PCA2 projection to the initial data
val ldaTexturesDFA = ldaModel.transform(whitener.transform(pcaModel
                             .transform(scaledTexturesDFA)))
                             .select("label","ldaFeatures")

Nous pouvons ensuite réaliser une visualisation avec Vegas :

import org.apache.spark.ml.linalg.Vector
import org.apache.spark.sql.functions.udf
// Extrait le premier scalaire d'un vecteur et le place dans une colonne
val first = udf((v: Vector) => v.toArray(0))
// Extrait le deuxième scalaire d'un vecteur et le place dans une colonne
val second = udf((v: Vector) => v.toArray(1))

val coordinates = ldaTexturesDFA.withColumn("x", first($"ldaFeatures")).withColumn("y", second($"ldaFeatures"))
coordinates.printSchema()

implicit val render = vegas.render.ShowHTML(s => print("%html " + s))

import vegas._
import vegas.data.External._
import vegas.sparkExt._
Vegas("AFD 2D").withDataFrame(coordinates).mark(Point).encodeX("x", Quant).encodeY("y", Quant).encodeColor("label", Nom).show

Les données appartenant à une même classe sont représentées avec une même couleur. Vu le nombre élevé de classes (11), certaines nuances sont proches entre classes différentes. Vous pouvez trouver une projection des mêmes données sur les 2 premiers axes discriminants dans la Fig. 13 de cette section de RCP216.