création svd et pca

This commit is contained in:
Pierre-Edouard Portier 2022-03-13 22:46:20 +01:00
parent cb4a0b3ef3
commit 183e0e0559
3 changed files with 186 additions and 4 deletions

3
.gitignore vendored
View File

@ -16,6 +16,9 @@
04_validation_croisee_slides.pdf 04_validation_croisee_slides.pdf
05_presentation_svd.pdf 05_presentation_svd.pdf
05_presentation_svd_slides.pdf 05_presentation_svd_slides.pdf
05_b_svd_pca.pdf
05_b_svd_pca_cache/
05_b_svd_pca_files/
06_matrice_definie_non_negative.pdf 06_matrice_definie_non_negative.pdf
07_multiplicateurs_de_lagrange.pdf 07_multiplicateurs_de_lagrange.pdf
08_derivation_svd.pdf 08_derivation_svd.pdf

179
05_b_svd_pca.Rmd Normal file
View File

@ -0,0 +1,179 @@
---
title: "05-b SVD et Analyse en Composantes Principales"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Projection sur les axes principaux
## SVD
Soit une matrice de données $\mathbf{X}\in\mathbb{R}^{n \times p}$ dont la décomposition en valeurs singlières est :
$$\mathbf{X} = \sum_{\alpha}\sqrt{\lambda_\alpha}\mathbf{u_\alpha}\mathbf{v_\alpha}^T \quad\equiv\quad \mathbf{X} = \mathbf{U} \mathbf{D} \mathbf{V}^T $$
## Facteurs
Nous avons montré que la projection orthogonale des lignes $\mathbf{x_i}$ de $\mathbf{X}$ sur $\mathbf{v_1},\dots,\mathbf{v_k}$ est la meilleure approximation $k$-dimensionnelle de $\mathbf{X}$ au sens de la minimisation des résidus au carré. Les axes de vecteurs directeurs $\mathbf{v_\alpha}$ sont appelés les \emph{axes principaux}. Les coordonnées des observations sur les axes principaux sont appelées \emph{facteurs}. Notons $\mathbf{F}$ la matrice dont les lignes sont les facteurs :
$$\mathbf{F} \triangleq \mathbf{X}\mathbf{V} = \mathbf{U} \mathbf{D} \mathbf{V}^T\mathbf{V} = \mathbf{U} \mathbf{D}$$
La covariance des facteurs est :
$$\left(\mathbf{U} \mathbf{D}\right)^T\left(\mathbf{U} \mathbf{D}\right) = \mathbf{D}^2$$
Nous vérifions ainsi que, par construction, les facteurs sont orthogonaux et que $\lambda_\alpha$ est la somme des facteurs au carré sur l'axe $\mathbf{v_\alpha}$, autrement dit, la variance expliquéee par cet axe.
## Contributions des observations
Nous pouvons mesurer la contribution $CTR_{i,\alpha}$ d'une observation $\mathbf{x_i}$ à la variance expliquée par $\mathbf{v_\alpha}$ :
$$CTR_{i,\alpha} = \frac{f_{i,\alpha}^2}{\sum_i f_{i,\alpha}^2} = \frac{f_{i,\alpha}^2}{\lambda_\alpha}$$
Puisque $\sum_i CTR_{i,\alpha} = 1$, nous pouvons considérer, de façon heuristique, que les observations qui contribuent le plus à expliquer l'axe $\mathbf{v_\alpha}$ sont telles que $CTR_{i,\alpha} \geq 1/n$. Les observations dont les contributions sont les plus importantes et de signes opposés peuvent permettre d'interpréter un axe en fonction de l'opposition de ses pôles.
## Contributions des axes
Nous pouvons similairement mesurer combien l'axe $\mathbf{v_\alpha}$ contribue à expliquer l'écart au centre de gravité d'une observation $\mathbf{x_i}$ :
$$COS2_{i,\alpha} = \frac{f_{i,\alpha}^2}{\sum_\alpha f_{i,\alpha}^2} = \frac{f_{i,\alpha}^2}{d_{i,g}^2}$$
Avec $d_{i,g}^2$ le carré de la distance de l'observation $\mathbf{x_i}$ au centre de gravité $g$. $d_{i,g}^2 = \sum_j \left(x_{i,j}-g_j\right)^2$. Si les données sont centrées alors $d_{i,g}^2 = \sum_j x_{i,j}^2$. La somme des carrés des distances au centre de gravité pour toutes les observations est égale à la variance totale des données, ou inertie totale : $\sum_i d_{i,j}^2 = \mathcal{I} = \sum_\alpha \lambda_\alpha$.
## Contribution des variables aux axes principaux
Les axes principaux $\mathbf{v_\alpha}$ sont des combinaisons linéaires des variables initiales. Lorsque la matrice initiale est centrée et réduite, les éléments de $(\lambda_\alpha/\sqrt{n-1})\mathbf{v_\alpha}$ représentent les corrélations entre les variables initiales et l'axe $\mathbf{v_\alpha}$. Ainsi, nous pouvons mesurer la contribution des variables initiales à l'expression de la variance expliquée par chaque axe principal (l'expression est au carré pour que la somme des contributions des variables à l'axe $\mathbf{v_\alpha}$ soit égale à $1$) :
$$VARCTR_{j,\alpha} = \left(\frac{\lambda_\alpha}{\sqrt{n-1}}\mathbf{v_\alpha}\right)^2$$
# Exemple
Considérons le jeu de données `abalone` introduit dans un précédent module. Nous retirons les observations pour lesquelles la variable `height` est nulle.
```{r cache=TRUE}
abalone.cols = c("sex", "length", "diameter", "height", "whole.wt",
"shucked.wt", "viscera.wt", "shell.wt", "rings")
url <- 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data'
abalone <- read.table(url, sep=",", row.names=NULL, col.names=abalone.cols,
nrows=4177)
abalone <- subset(abalone, height!=0)
```
Nous sélectionnons les variables explicatives numériques dans une matrice $\mathbf{X}$.
```{r}
X <- abalone[,c("length", "diameter", "height", "whole.wt","shucked.wt", "viscera.wt",
"shell.wt")]
```
Nous écrivons une fonction pour standardiser $\mathbf{X}$.
```{r}
standard <- function(X) {
meanX <- apply(X, 2, mean);
sdX <- apply(X, 2, sd);
Y <- sweep(X, 2, meanX, FUN = "-");
Y <- sweep(Y, 2, sdX, FUN = "/");
res <- list("X" = Y, "meanX" = meanX, "sdX" = sdX)
}
```
```{r warning=FALSE}
set.seed(1123)
std.X <- standard(X)
A <- std.X$X
nbClust <- 100
clst <- kmeans(A, nbClust, nstart = 25)
```
Avec le code ci-dessus, nous utilisons l'algorithme `k-means` pour découvrir `r nbClust` clusters à partir des données standardisées.
Nous appliquons une décomposition en valeurs singulières sur les centres de clusters standardisés.
```{r}
std.centers <- standard(clst$centers)
C <- std.centers$X
svd.C <- svd(C)
```
Nous calculons le pourcentage de variance expliquée par chaque axe principal.
```{r}
round((svd.C$d^2 / sum(svd.C$d^2))*100, 2)
```
Nous calculons les facteurs et les contributions des observations aux axes principaux.
```{r}
fact <- svd.C$u %*% diag(svd.C$d)
fact2 <- fact^2
ctr <- sweep(fact2, 2, colSums(fact2), "/")
```
Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes.
```{r}
n <- dim(fact)[1]
d1BestCtr <- which(ctr[,1] > 1/n)
d2BestCtr <- which(ctr[,2] > 1/n)
d1d2BestCtr <- union(d1BestCtr,d2BestCtr)
plot(fact[d1d2BestCtr,1], fact[d1d2BestCtr,2], pch="")
text(fact[d1d2BestCtr,1], fact[d1d2BestCtr,2], d1d2BestCtr, cex=0.8)
```
```{r}
f2.max <- which.max(abs(fact[,2]))
f2.max.size <- clst$size[f2.max]
f2.max.name <- names(clst$cluster[clst$cluster==f2.max])
```
Le cluster qui explique le plus la variance du deuxième axe principal semble anormalement important. Il contribue à expliquer `r round(ctr[f2.max,2]*100, 2)` % de la variance du second axe. Il est composé de seulement `r f2.max.size` élément, en retournant au jeu de données original :
```{r}
abalone[f2.max.name,]
```
Nous retrouvons l'observation anormale déjà identifiée dans une précédente analyse. Nous recommençons l'analyse en le retirant :
```{r, warning=FALSE}
abalone <- abalone[!(row.names(abalone) %in% f2.max.name),]
X <- abalone[,c("length", "diameter", "height", "whole.wt","shucked.wt", "viscera.wt",
"shell.wt")]
std.X <- standard(X)
A <- std.X$X
nbClust <- 100
clst <- kmeans(A, nbClust, nstart = 25)
std.centers <- standard(clst$centers)
C <- std.centers$X
svd.C <- svd(C)
round((svd.C$d^2 / sum(svd.C$d^2))*100, 2)
```
Après cette correction, nous voyons qu'une part encore plus importante de la variance est expliquée par le premier axe principal. Cela peut nous indiquer que les variables explicatives sont très corrélées.
Nous calculons à nouveau les facteurs et les contributions des observations aux axes principaux. Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes.
```{r}
fact <- svd.C$u %*% diag(svd.C$d)
fact2 <- fact^2
ctr <- sweep(fact2, 2, colSums(fact2), "/")
n <- dim(fact)[1]
d1BestCtr <- which(ctr[,1] > 1/n)
d2BestCtr <- which(ctr[,2] > 1/n)
d1d2BestCtr <- union(d1BestCtr,d2BestCtr)
plot(fact[d1d2BestCtr,1], fact[d1d2BestCtr,2], pch="")
text(fact[d1d2BestCtr,1], fact[d1d2BestCtr,2], d1d2BestCtr, cex=0.8)
```
Le cluster `43` `r intrig<-43` est intrigant. Sa taille est : `r clst$size[intrig]`. C'est encore un singleton qui correspond à l'observation `r names(clst$cluster[clst$cluster==intrig])` du jeu de données original. Nous découvrons qu'il s'agit sans doute d'une anomalie pour la variable `height` :
```{r}
intrig.name <- names(clst$cluster[clst$cluster==intrig])
abalone[intrig.name,]
```
```{r}
boxplot(abalone$height)
```
Profitons-en pour calculer les contributions des axes principaux aux écarts au centre d'inertie des observations.
```{r}
cos2 <- sweep(fact2, 1, rowSums(fact2), "/")
round(cos2[intrig,]*100,2)
```
Nous vérifions que cette observation, sans doute anormale, est presque entièrement expliquée par les deux premiers axes principaux. Nous ne prenons donc pas de risque à considérer comme significatif son écart à l'origine dans le plan formé par ces deux axes.
Observons aussi les contributions des variables aux axes principaux :
```{r}
varctr <- ((svd.C$v %*% diag(svd.C$d))/sqrt(n-1))^2
rownames(varctr) <- colnames(X)
round(varctr,2)
```
Nous voyons à nouveau que toutes les variables sont très corrélées car elles contribuent toutes fortement à la variance expliquée par le premier axe principal.

View File

@ -11,9 +11,9 @@ source("05_presentation_svd.R", local = knitr::knit_global())
# Géométrie d'un tableau de données # Géométrie d'un tableau de données
- $\mathbf{X} \in \mathcal{R}^{n \times p}$ ; $x_{ij} \in \mathcal{R}$ - $\mathbf{X} \in \mathbb{R}^{n \times p}$ ; $x_{ij} \in \mathbb{R}$
- Nuage de $n$ points $\mathbf{x_i}$ dans $\mathcal{R}^p$ - Nuage de $n$ points $\mathbf{x_i}$ dans $\mathbb{R}^p$
- Nuage de $p$ points $\mathbf{x_j}$ dans $\mathcal{R}^n$ - Nuage de $p$ points $\mathbf{x_j}$ dans $\mathbb{R}^n$
- Projeter le nuage des points $\mathbf{x_i}$ sur $\mathcal{H} \subset \mathbb{R}^P$ - Projeter le nuage des points $\mathbf{x_i}$ sur $\mathcal{H} \subset \mathbb{R}^P$
- Minimiser les déformations - Minimiser les déformations
@ -37,7 +37,7 @@ knitr::include_graphics("images/svd1.jpeg")
- $\mathbf{X}^T \mathbf{X} \mathbf{v_1} = \lambda_1 \mathbf{v_1}$ avec $\lambda_1$ la plus grande valeur propre - $\mathbf{X}^T \mathbf{X} \mathbf{v_1} = \lambda_1 \mathbf{v_1}$ avec $\lambda_1$ la plus grande valeur propre
- $\mathbf{X}^T \mathbf{X} \mathbf{v_2} = \lambda_2 \mathbf{v_2}$ avec $\lambda_2$ la seconde plus grande valeur propre - $\mathbf{X}^T \mathbf{X} \mathbf{v_2} = \lambda_2 \mathbf{v_2}$ avec $\lambda_2$ la seconde plus grande valeur propre
- $\mathbf{v_1}, \mathbf{v_2},\dots \mathbf{v_k}$ forment une base de $\mathcal{H}_k$ - $\mathbf{v_1}, \mathbf{v_2},\dots \mathbf{v_k}$ forment une base de $\mathcal{H}_k$
- La projection du nuage des $n$ points $\mathbf{x_i}$ sur $\mathcal{H}_k$ minimise les carrés des écarts à $\mathcal{R}^p$ - La projection du nuage des $n$ points $\mathbf{x_i}$ sur $\mathcal{H}_k$ minimise les carrés des écarts à $\mathbb{R}^p$
- Symétriquement pour les vecteurs $\mathbf{u_1}, \mathbf{u_2},\dots \mathbf{u_k}$, etc. - Symétriquement pour les vecteurs $\mathbf{u_1}, \mathbf{u_2},\dots \mathbf{u_k}$, etc.
# Décomposition en valeurs singulières # Décomposition en valeurs singulières