Décomposition en valeurs singulières SVD d'une matrice rectangulaire comme introduction à l'Analyse en Composantes Principales (ACP); illustration sur des données fictives avec R puis avec le package FactoMineR; application à des données élémentaires: courbes annuelles de températures moyennes mensuelles de 36 villes françaises.
Les différents travaux et analyses proposés tout au long de ces documents sont largement explicités. Les commandes en R ou Python toutes fournies. L'important n'est pas de trouver la bonne syntaxe des commandes ni de finir au plus vite mais de réfléchir sur les méthodes, leurs conditions d'applications, les résultats obtenus. L'apprentissage de ces logiciels et de leur programmation est dispensés dans d'autres tutoriels.
Répondre aux questions Q tout au long de la réalisation.
Cette section se propose d'illustrer les principaux résultats d'algèbre linéaire utiles en exploration statistique multidimensionnelle. Ceci concerne les valeurs et vecteurs propres de matrices symétriques en lien avec la décomposition en valeur singulière ou SVD d'une matrice rectangulaire $n\times p$ pour en faire une approximation par une matrice de mêmes dimensions mais de rang inférieur.
Décomposition en valeurs singulières d'une matrice rectangulaire relativement à des métriques classiques définies par la matrice identité.
Q Ecrire la formule de la décomposition en valeurs singulières d'une matrice rectangulaire.
Q De quelles matrice les vecteurs singuliers sont vecteurs propres?
In [ ]:
# Génération d'un matrice n x p
X=matrix(runif(100),20,5)
# SVD
res=svd(X)
# Valeurs singulières
res$d
In [ ]:
# Vérifier l'orthonormalité des vecteurs
t(res$u)%*%res$u
In [ ]:
t(res$v)%*%res$v
In [ ]:
# Vérifier la reconstruction de X
# à l'erreur machine près
X-res$u%*%diag(res$d)%*%t(res$v)
Q Comparer ci-dessous les valeurs propres des matrices X'X et XX'. Que dire des dimensions, du rang de ces matrices, de la multiplicité de la valeur propre nulle?
Q Comparer avec les valeurs propres et les valeurs singulières de X.
In [ ]:
# Valeurs et vecteurs propres
dec1=eigen(t(X)%*%X)
dec2=eigen(X%*%t(X))
dec1$values
In [ ]:
dec2$values
In [ ]:
sqrt(dec1$values)
In [ ]:
U=dec2$vectors
V=dec1$vectors
# Orthonormalité des vecteurs
t(V)%*%V
In [ ]:
t(U)%*%U
Q Vérifier la bonne cohérence des dimensions de ces vecteurs puis comparer les vecteurs singuliers à droite et à gauche de X avec les vecteurs propres de ces matrices. D'où viennent les différences ?
In [ ]:
V-res$v
In [ ]:
U[,1:5]-res$u
Q Vérifier que les premiers termes de la SVD sont la meilleure approximation de X par une matrice de rang inférieure.
In [ ]:
# approximation de rang 4
Xhat=res$u[,1:4]%*%diag(res$d[1:4])%*%t(res$v[,1:4])
# calcul de la norme ||X-Xhat||^2
sum((X-Xhat)**2)
Q Comparer cette norme avec les carrés des valeurs singulières, que vaudrait cette norme avec une approximation de rang 3?
Q Retrouver ci-dessous la matrice des vecteurs singuliers à droite à partir de celle des vecteurs singuliers à gauche et réciproquement.
In [ ]:
# Vecteurs singuliers à gauche:
Ug=res$u
# vecteurs singuliers à droites
# calculés à partir de U
Vd=t(X)%*%Ug
# Vérifier l'orthogonalité
t(Vd)%*%Vd
mais pas l'orthonormalité
Q Que vaut la diagonale ?
In [ ]:
diag(t(Vd)%*%Vd)-res$d**2
Il est nécessaire de normer les vecteurs.
In [ ]:
Vd=Vd%*%diag(1/res$d)
t(Vd)%*%Vd
In [ ]:
# Vérifier
Vd-res$v
Q Retrouver de façon analogue les vecteurs singuliers Ug à partir de ceux Vd en prémultipliant par la matrice X.
Les calculs précédents concernent des matrices (applications linéaires) dans des espaces euclidiens munis de métriques usuelles définies par les matrices identités. Pour calculer une SVD généralisée par rapport à des métriques quelconques, celà nécessite d'introduire des changements de métriques car la SVCD généralisée n'est pas prévue dans les logiciels.
Comme en ACP, on considère une métrique dans l'espace vectoriel $\boldsymbol{R}^p (p=5)$ définie par une matrice $\boldsymbol{M}$ symétrique définie positive et une métrique de matrice diagonale $\boldsymbol{D}$ sur l'expace vectoriel $\boldsymbol{R}^n (n=20)$. La SVD de $\boldsymbol{X}$ relativement à $\boldsymbol{M}$ et $\boldsymbol{D}$ est calculée à partir de celle classique de $\boldsymbol{D}^{1/2}\boldsymbol{X}\boldsymbol{M}^{1/2}$ relativement aux métriques de matrice identité; $\boldsymbol{M}^{1/2}$ est la racine carrée de $\boldsymbol{M}$.
En effet, si $\boldsymbol{D}^{1/2}\boldsymbol{X}\boldsymbol{M}^{1/2}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}' \quad \text{avec} \quad \boldsymbol{U'U}=\boldsymbol{I} \quad \text{et} \quad \boldsymbol{V'V}=\boldsymbol{I}$
alors en prémultipliant à gauche et multipliant à droite:
$ \boldsymbol{X}=\boldsymbol{D}^{-1/2}\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}'\boldsymbol{M}^{-1/2} = \left(\boldsymbol{D}^{-1/2}\boldsymbol{U}\right)\boldsymbol{\Lambda}\left(\boldsymbol{M}^{-1/2}\boldsymbol{V}\right)' $
et $\boldsymbol{D}^{-1/2}\boldsymbol{U}$ définit les vecteurs $\boldsymbol{D}$-orthonormés singuliers à gauche et $\boldsymbol{M}^{-1/2}\boldsymbol{V}$ ceux à droite de la SVD$(\boldsymbol{X,M,D})$ car ils vérifient bien: $ \left(\boldsymbol{D}^{-1/2}\boldsymbol{U}\right)'\boldsymbol{D}\left(\boldsymbol{D}^{-1/2}\boldsymbol{U}\right)=\boldsymbol{I} \quad\text{et}\quad \left(\boldsymbol{M}^{-1/2}\boldsymbol{V}\right)'\boldsymbol{M}\left(\boldsymbol{M}^{-1/2}\boldsymbol{V}\right)=\boldsymbol{I}. $
Construction des différentes matrices.
In [ ]:
# Matrice M symétrique définie positive
H=matrix(rnorm(25),5,5)
M=t(H)%*%H
# Décompostion de M
s=eigen(M);l=s$values;v=s$vectors
l; v
In [ ]:
# Racine de M
Mr=v%*%diag(sqrt(l))%*%t(v)
# Inverse de racine de M
Mi=v%*%diag(1/sqrt(l))%*%t(v)
# vérifications
Mr%*%Mi;M-Mr%*%Mr
In [ ]:
# Construction de D
w=runif(20); w=w/sum(w)
D=diag(w)
Dr=diag(sqrt(w)); Di=diag(1/sqrt(w))
Construction de la SVD$(\boldsymbol{X,M,D})$ et vérifications.
In [ ]:
# SVD de l'image de X (changement de métrique)
s=svd(Dr%*%X%*%Mr)
U=s$u;L=s$d;V=s$v
# Nouveaux vecteurs singuliers
Ud=Di%*%U
Vm=Mi%*%V
# Vérifications
t(Vm)%*%M%*%Vm
In [ ]:
t(Ud)%*%D%*%Ud
X-Ud%*%diag(L)%*%t(Vm)
L'objectif est de bien comprendre quels sont les résultats / matrices produites par les différentes façons d'aborder l'analyse en composantes principales.
Les données sont celles de l'exemple introduction à l'ACP: les notes en maths, français, physique et anglais de 9 lycéens virtuels.
In [ ]:
# Matrice des données
note=matrix(c(6,6,5,5.5,8,8,8,8,
6,7,11,9.5,14.5,14.5,15.5,15,
14,14,12,12.5,11,10,5.5,7,
5.5,7,14,11.5,13,12.5,8.5,9.5,
9,9.5,12.5,12),nrow=9,byrow=TRUE)
note=data.frame(note)
nomi=c("jean","alai","anni","moni",
"didi","andr","pier","brig","evel")
nomv=c("Math","Phys","Fran","Angl")
dimnames(note)[[1]]=nomi
dimnames(note)[[2]]=nomv
# Vérification
note
Statistiques élémentaires :
In [ ]:
summary(note)
options(repr.plot.width=5, repr.plot.height=3)
boxplot(note)
In [ ]:
X=as.matrix(note)
# Matrice des variances covariances
var(X)
In [ ]:
# Matrice des corrélations
cor(X)
Q Commenter la matrice des corrélations ci-dessus.
Q Quel est le graphique ci-dessous? Comment l'interpréter?
In [ ]:
pairs(X,pch=19, col="blue")
Tout langage matricielle permet de construire très facilement une ACP. Attention, la formule de variance utilise $(n-1)$ (estimateur sans biais) au lieu de $n$ comme diviseur afin de retrouver les résultats des fonctions R. Un logiciel français peut donc fournir des résultats légèrement différents. Un des objectifs des cette approche élémentaire sur des données jouet est justement de faire le tri.
In [ ]:
# Matrice centrée
Xb=scale(X,scale=F)
Xb
In [ ]:
# Covariance
S=t(Xb)%*%Xb/8
# SVD
res=svd(Xb)
# Matrices des vecteurs propres
U=res$u;V=res$v
# Valeurs singulières
vs=res$d
Valeurs singulières de $(\bar{X}, I_p, I_n)$ puis avec $1/n I_n$ sivi de $1/(n-1) I_n$
In [ ]:
vs
In [ ]:
vs/sqrt(8)
In [ ]:
vs/sqrt(9)
Valeurs propres avec $n-1$ comme diviseur
In [ ]:
(vs**2)/8
Valeurs propres avec $n$ comme diviseur
In [ ]:
(vs**2)/9
Vecteurs propres de $\bar{X}'\bar{X}/(n-1)$ ou $/n$ ou vecteurs singuliers à droite de $\bar{X}$
In [ ]:
V
In [ ]:
# Valeurs propres
L=res$d*res$d/8; pct=L/sum(L)
L;pct
Q Quelle signification statistique donner aux valeurs propres?
In [ ]:
# Composantes principales
C=Xb%*%V
Construction de représentations graphiques rudimentaires.
Q Quel est le raphe ci-dessous? Quelle conclusion en tirer?
In [ ]:
boxplot(as.data.frame(C))
In [ ]:
# Coordonnées "isométriques lignes"
options(repr.plot.width=5, repr.plot.height=4)
plot(C,type="n")
text(C,nomi)
abline(h=0,v=0)
Q Quelle insuffisance présente ce graphique? Quel choix d'options tente de la corriger?
In [ ]:
# Coordonnées "isométriques colonnes"
plot(V[,1]*sqrt(L[1]),V[,2]*sqrt(L[2]),type="n")
text(V[,1]*sqrt(L[1]),V[,2]*sqrt(L[2]),nomv)
abline(h=0,v=0)
Plus d'effort est nécessaire pour produire des graphiques lisibles. Autant utiliser les fonctions R adaptées. D'abord avec la fonction princomp.
Attention: L'utilisation de fonctions prédéfinies simplifie le travail mais nécessite de bien en contrôler les options part défaut afin de savoir y lire les bons résultats. Deux fonctions sont disponibles dans R: princomp et prcomp.
princomp
In [ ]:
res1.acp=princomp(note)
summary(res1.acp)
In [ ]:
attributes(res1.acp)
In [ ]:
res1.acp$sdev
Q Quel est le diviseur de la variance pour princomp.
Vecteurs propres au signe près.
In [ ]:
res1.acp$loadings
In [ ]:
plot(res1.acp)
In [ ]:
biplot(res1.acp)
In [ ]:
res2.acp=prcomp(note)
summary(res2.acp)
In [ ]:
attributes(res2.acp)
In [ ]:
res2.acp$sdev
Q Quel est le diviseur de la variance pour prcomp.
Vecteurs propres
In [ ]:
res2.acp$rotation
In [ ]:
plot(res2.acp)
In [ ]:
biplot(res2.acp)
Q Comparer les aides de ces fonctions pour comprendre d'où viennent les différences. Attention pour princomp, le nombre de lignes ($n$) doit être plus grand que le nombre de colonnes $p$; ce n'est pas une contrainte pour prcomp.
Q Que sont les différents résultats de ces deux fonctions : sdev rotation center scale x de prcomp et sdev loadings center scores de princomp par rapport aux matrices U, V, L, C précédentes.
Ce package donne accès à la plupart des méthodes factorielles et de classification non supervisée multidimensionnelles. Son utilisation nécessite, si ce n'est déjà fait, une installation préalable par la commande install.packages("FactoMineR") ou par l'utilisation des menus de la fenêtre inférieure droite de RStudio.
In [ ]:
library(FactoMineR)
acp=PCA(note, graph=FALSE)
Ce package fournit de très (un peu trop) nombreux résultats et des graphiques dont le flux n'est pas pas compatible avec jupyter. Il est préférable de commander chaque résultat ou graphique dans cet environnement.
In [ ]:
print(acp)
In [ ]:
acp$eig
Q Acp est-elle réduite par défaut?
In [ ]:
acp=PCA(note, graph=FALSE, scale.unit=FALSE)
acp$eig
In [ ]:
acp$svd$vs
Q Quel est le diviseur de la variance?
Q Comment retrouver la matrice $V$ des vecteus propres à partir des coordonnées des variables.
In [ ]:
acp$var$coord%*%diag(1/acp$svd$vs)
In [ ]:
plot(acp,title="")
In [ ]:
plot(acp, choix="var",title="")
In [ ]:
acp=PCA(note, graph=FALSE)
plot(acp,choix="var")
Q Quel est le cercle qui apparaît sur ce graphique. Quelle interprétation en tirer?
En résumé: le diviseur de la variance dans princomp et PCA de FactoMineR est $n$, celui de prcomp est $n-1$. La fonction PCA de FactoMineR réduit par défaut mais pas princomp ni prcomp.
Les données étudiées sont celles du fichier tempR.dat disponibles dans le même dépôt que ce calepin. Il contient les moyennes, entre 1931 et 1960, des températures mensuelles moyennes de 36 villes françaises. La première variable correspond au nom de la ville (4 caractères), les 12 suivantes représentent chacune un mois de l'année (source : Mémorial de la Météorologie nationale). Une moyenne journalière est la moyenne des températures min et max permettant ensuite de calculer la moyenne mensuelle. Moyenner ensuite sur 10 ces valeurs conduit à des courbes relativement régulières.
In [ ]:
# lire les données dans R
temp=read.table("Data/tempR.dat")
# vérification et stat élémentaires
summary(temp)
In [ ]:
# Ensemble des courbes
options(repr.plot.width=4, repr.plot.height=4)
ts.plot(t(temp),col="blue")
In [ ]:
# toutes les distributions
boxplot(temp)
Q Que dire des dispersions, de l'homogénéité des unités et des variances?
Q Conséquence pour l'ACP?
In [ ]:
# Corrélations
cor(temp)
In [ ]:
plot(temp$janv,temp$fevr)
In [ ]:
plot(temp$janv,temp$juin)
In [ ]:
options(repr.plot.width=8, repr.plot.height=8)
pairs(temp)
In [ ]:
acp=princomp(temp)
summary(acp)
options(repr.plot.width=4, repr.plot.height=3)
plot(acp)
In [ ]:
boxplot(data.frame(acp$scores))
In [ ]:
options(repr.plot.width=6, repr.plot.height=6)
biplot(acp)
Q Analyser les échelles des axes.
Q Expliciter le choix entre ACP réduite ou non, comparer les différences.
Q Combien d'axes faut-il retenir ? Justifier.
Q Identifier la ville atypique de l'axe 2. Que faire ?
Q Interprétation des axes.
Q Commenter la position d'Embrun sur les graphiques.
Les données sont ici très particulières, des courbes fonction du temps. En conséquence, les vecteurs propres le sont également et les courbes sont décomposées sur cette base de fonctions discrétisées.
Vecteurs ou plutôt fonctions propres de l'ACP.
In [ ]:
options(repr.plot.width=5, repr.plot.height=5)
plot.ts(acp$loadings[,1:6], main="Fonctions propres")
Cette librairie apporte des compléments intéressants (qualité et options des graphiques, gestion des variables manquantes) et surtout elle vient particulièrement compléter les fonctions de base de R pour l'analyse des variables qualitatives. Voici les principaux résultats de l'ACP.
Q Comparer avec les résultats numériques précédemment obtenus.
In [ ]:
library(FactoMineR)
acp=PCA(temp, scale.unit=FALSE,ncp=12,graph=F)
options(repr.plot.width=4, repr.plot.height=3)
barplot(acp$eig[,1])
In [ ]:
boxplot(acp$ind$coord)
In [ ]:
acp$svd$V
In [ ]:
dimdesc(acp,axes=c(1,2))
Attention Cette librairie ajoute à ces techniques exploratoires des résultats d'inférence statistique: p-valeurs de test, ellipse de confiance... supposant implicitement un modèle probabiliste: distributions gaussiennes multidimensionnelles; ils sont à utiliser avec grande prudence, plus comme des indicateurs que comme des aides formelles à la décision.
In [ ]:
options(repr.plot.width=5, repr.plot.height=5)
plot(acp)
In [ ]:
plot(acp,choix="var")
In [ ]:
acp=PCA(temp, scale.unit=TRUE,ncp=12,graph=F)
plot(acp,choix="var")
Q Quelle différence entre les deux graphiques? Quelle différence entre les deux ACP?
Q Pourquoi les deux réprésentations sont-elles finalement très similaires?
Il s'agit également de résoudre la question concernant l'observation atypique sur l'axe 2. Faut-il la conserver? Cette question est abordée en considérant deux ACPs, celle avec et sans ce point afin de s'assurer que sa suppression ne perturbe pas trop les premiers axes, notamment le 2ème.
Q Comment l'observation atypique et-elle exclue?
In [ ]:
acp=PCA(temp, scale.unit=TRUE,ncp=12,graph=F,ind.sup=8)
plot(acp,choix="var")
Q L'interprétation des axes a-t-elle été modifiée? Que faire de l'observation atypique?