p6e7p7/intro_to_ml
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modify the individual chapters to create a bookdown book

Browse Source
This commit is contained in:
Pierre-Edouard Portier 2022-04-16 11:52:56 +02:00
parent 2f06eaed4b
commit 60ecf4ee1c
31 changed files with 411 additions and 470 deletions

8
01_intro.R

@ -11,9 +11,9 @@ gendat <- function(n, sd) {
# sd: écart type d'un bruit gaussien de moyenne nulle
X = runif(n)
Y = f(X) + rnorm(n, mean=0, sd=sd)
dim(X) <- c(n,1) # en général chaque observation est décrite par plusieurs
# variables et X est une matrice avec autant de lignes que
# d'observations et autant de colonnes que de variables. Sur
dim(X) <- c(n,1) # en général chaque observation est décrite par plusieurs
# variables et X est une matrice avec autant de lignes que
# d'observations et autant de colonnes que de variables. Sur
# notre exemple, chaque observation n'est décrite que par une
# seule variable.
list(X = X, Y = Y)
@ -47,4 +47,4 @@ polyeval <- function(coef,x) {
pltpoly <- function(coef) {
xs = seq(0,1,length.out=100)
lines(xs, polyeval(coef,xs), lty="dotted")
}
}

32
01_intro.Rmd

@ -1,17 +1,10 @@
---
title: "01 Introduction"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Introduction
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
```
# Contexte
## Contexte
A l’occasion des trois premières révolutions industrielles, des tâches, auparavant réservées au travail manuel de l’Homme, ont été automatisées. Il semble envisageable d’associer au tournant du 21ème siècle une quatrième révolution portée par l’automatisation de la capacité à prédire, essentielle pour le processus de prise de décision dans les secteurs de l’industrie, du commerce et des services.
@ -21,7 +14,7 @@ Par ailleurs, cette automatisation des prédictions s’accompagne d’un renouv
Ainsi, ce cours introductif à l’apprentissage automatique a pour objectif d’offrir des connaissances fondamentales et des compétences pratiques qui aideront l'ingénieur à tenir ce rôle essentiel.
# Objectifs
## Objectifs
La discipline de l’apprentissage automatique, ou machine learning, élabore des algorithmes et des méthodes pour découvrir des régularités dans des données multidimensionnelles afin, entre autres, d’automatiser la prédiction. Elle peut se subdiviser en trois catégories. D’abord, l’apprentissage non (ou semi) supervisé qui s’attache à découvrir des structures dans les données non étiquetées à travers des approches comme le clustering, la réduction dimensionnelle, les modèles génératifs… Ensuite, l’apprentissage par renforcement, dans le cadre duquel un agent interagit avec son environnement en adaptant son comportement pour maximiser une fonction de récompense. Enfin, l’apprentissage supervisé, qui fait l’objet de ce module, a quant à lui pour objectif d’apprendre à prédire l’association entre un objet décrit selon plusieurs dimensions et une étiquette.
@ -31,13 +24,13 @@ Ainsi, ce cours a pour objectif d’introduire quelques concepts fondamentaux de
Pour faciliter l’acquisition des connaissances, le cours est accompagné de notebooks manipulables, rédigés dans le langage de programmation R. Ces mises en pratique systématiques doivent permettre de faire le lien entre des concepts fondamentaux et leur application dans des projets d’analyse de données.
# Ajustement de courbe
## Ajustement de courbe
$\mathbf{x}^{(1)} \dots \mathbf{x}^{(n)}$ sont des vecteurs de $\mathbb{R}^p$ associés aux valeurs, aussi appelées étiquettes, $y^{(1)} \dots y^{(n)}$ de $\mathbb{R}$. Nous cherchons une fonction $f(\mathbf{x}) : \mathbb{R}^p \rightarrow \mathbb{R}$ qui modélise la relation entre les observations $\mathbf{x}$ et les étiquettes $y$.
La fonction $f$ peut avoir une forme paramétrique, comme par exemple :
\[
f(\mathbf{x}) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p
f(\mathbf{x}) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p
\]
Si $n=p+1$, les paramètres $\beta_0, \beta_1, \dots \beta_p$ sont solution d'un système linéaire :
\[
@ -61,7 +54,7 @@ Ce système s'écrit également sous forme matricielle :
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
\end{array} \right)
\]
Chaque ligne $i$ de la matrice du terme de gauche de l'égalité ci-dessus est le vecteur ligne $\mathbf{x}^{(i)T}$ avec l'addition d'un premier terme constant qui correspond au paramètre $\beta_0$. En nommant cette matrice $\mathbf{X}^T$, le système linéaire ci-dessus s'écrit :
@ -82,17 +75,17 @@ Avec $n=p+1$ observations et les étiquettes associées $\left( x^{(k)}, y^{(k)}
1 & x^{(n)} & (x^{(n)})^2 & \dots & (x^{(n)})^p
\end{array} \right)
\left( \begin{array}{c}
\beta_0 \\ \beta_1 \\ \dots \\ \beta_p
\beta_0 \\ \beta_1 \\ \dots \\ \beta_p
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
\end{array} \right)
\]
La matrice du terme de gauche de l'égalité ci-dessus est traditionnellement appelée "matrice de Vandermonde".
# Interpolation polynomiale sur un jeu de données synthétique
## Interpolation polynomiale sur un jeu de données synthétique
Soit un exemple de fonction non-linéaire, $f(x) = e^x \times cos(2 \pi \times sin(\pi x))$, utilisée pour générer un jeu de données synthétique.
@ -115,4 +108,9 @@ plt(data,f)
pltpoly(coef)
```
Il est improbable que ce polynôme, passant exactement par chaque observation, puisse offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, nous cherchons à découvrir un polynôme de degré plus faible. Il ne passera pas exactement par toutes les observations, mais il prédira probablement mieux les étiquettes associées à de nouvelles observations.
Il est improbable que ce polynôme, passant exactement par chaque observation, puisse offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, nous cherchons à découvrir un polynôme de degré plus faible. Il ne passera pas exactement par toutes les observations, mais il prédira probablement mieux les étiquettes associées à de nouvelles observations.
## Annexe code source
```{r, code=readLines("01_intro.R"), eval=FALSE}
```

22
02_a_application_abalone.Rmd

@ -1,15 +1,6 @@
---
title: "02-a Application au jeu de données `abalone`"
output:
bookdown::pdf_document2:
number_section: yes
includes:
in_header: preamble.tex
toc: false
classoption: fleqn
---
# Application au jeu de données `abalone`
# Récupération du jeu de données
## Récupération du jeu de données
Téléchargeons le jeu de données [`abalone`](http://archive.ics.uci.edu/ml/datasets/Abalone) depuis le répertoire de l'Université de Californie Irvine (University of California Irvine, UCI, machine learning repository).
@ -38,7 +29,7 @@ abalone <- read.table(url, sep=",", row.names=NULL, col.names=abalone.cols,
str(abalone)
```
# Préparation du jeu de données
## Préparation du jeu de données
En cherchant des valeurs aberrantes, nous remarquons que, pour deux observations, la variable `height` a une valeur nulle.
```{r}
@ -68,7 +59,7 @@ vif(lm(rings ~ ., data = abalone))
Parmi les variables corrélées, nous proposons de conserver celles avec les scores VIF les plus faibles (c'est-à-dire celles dont la variance s'explique le moins par la variance d'autres variables).
# Modèle linéaire
## Modèle linéaire
Nous calculons un modèle linéaire.
```{r}
@ -80,12 +71,13 @@ Observons les résultats statistiques proposés par défaut.
summary(abalone_lm)
```
# Analyse graphique du modèle linéaire
## Analyse graphique du modèle linéaire
Affichons les résidus pour chaque observation.
```{r}
plot(abalone_lm, which = 1, id.n = 5)
```
L'observation d'identifiant $2052$ semble anormale. Nous observons que sa valeur `height` est anormalement élevée : `r abalone["2052","height"]`.
```{r}
abalone <- abalone[!(row.names(abalone) %in% c("2052")),]
@ -118,7 +110,7 @@ Observons les résidus en fonction des effets levier.
plot(abalone_lm, which = 5, id.n = 5)
```
# Interprétation du modèle linéaire
## Interprétation du modèle linéaire
Une augmentation de `length` de $1mm$ augmente la prédiction de $1.1$ anneaux avec un intervalle de confiance à $95\%$ qui vaut $[0.95,1.27]$.
De même pour les autres variables continues.

35
02_b_breviaire_proba_stat.Rmd

@ -1,17 +1,8 @@
---
title: "02-b Bréviaire proba/stat"
output:
bookdown::pdf_document2:
number_section: yes
includes:
in_header: preamble.tex
toc: false
classoption: fleqn
---
# Bréviaire proba/stat
# Probabilités
## Probabilités
## Espérance
### Espérance
$$
\begin{aligned}
@ -24,7 +15,7 @@ E[UV] &= E[U] E[V] \quad \text{si $U$ et $V$ sont indépendantes} \\
\end{aligned}
$$
## Variance
### Variance
$$
\begin{aligned}
@ -42,7 +33,7 @@ Var[aU+bV] &= a^2Var[U] + b^2Var[V] + 2 ab Cov[U,V] \\
\end{aligned}
$$
# Distribution normale
## Distribution normale
$$
\begin{aligned}
@ -53,7 +44,7 @@ X &\sim \mathcal{N}(\mu,\sigma^2) \quad \text{$X$ est distribué selon $\mathcal
\end{aligned}
$$
# Statistiques
## Statistiques
Pour des $X_i$ indépendantes identiquement distribuées, de moyenne $\mu$ et de variance $\sigma^2$ :
$$
@ -67,7 +58,7 @@ E[s^2] &= \frac{n-1}{n} \sigma^2 \quad \text{Estimateur biaisé de la variance.}
\end{aligned}
$$
# Théorème Central Limite (CLT)
## Théorème Central Limite (CLT)
Si $X_1,\dots,X_n$ sont des variables aléatoires indépendantes qui suivent la même distribution de moyenne $m$ et de variance $v^2$, alors $T=X_1+\dots+X_n$ suit approximativement $\mathcal{N}\left(nm,nv^2\right)$. L'approximation est en pratique bonne à partir d'environ $n=25$.
@ -84,7 +75,7 @@ $$P\left(| \bar{X} - \mu | \geq 3\sigma/\sqrt{n}\right) \leq \frac{1}{3^2} $$
Soit :
$$P\left(| \bar{X} - \mu | < 3\sigma/\sqrt{n}\right) \geq \frac{8}{9} = 0.8889 $$
# Intervalles de confiance
## Intervalles de confiance
C'est sur la base du CLT que nous déterminons des intervalles de confiance.
$$
@ -99,14 +90,14 @@ $$
Dans cette formule, il est juste de remplacer l'écart-type $\sigma$, inconnu, par son estimation $s$.
# Tests d'hypothèses
## Tests d'hypothèses
- Faire une hypothèse : $H_0 : \theta = c$
- Calculer $Z = \frac{\hat{\theta} - c}{\text{s.e.}(\hat{\theta})}$
- Rejeter $H_0$ à un niveau de risque $\alpha = 0.005$ si $|Z|>=1.96$
- La p-valeur est le plus faible niveau de risque auquel l'hypothèse serait rejetée.
# Application aux coefficients d'un modèle linéaire
## Application aux coefficients d'un modèle linéaire
Sous hypothèse d'un modèle génératif linéaire :
$$y_i = \mathbf{X}\boldsymbol\beta + \epsilon_i \quad , \quad \epsilon_i \sim \mathcal{N}(0,\sigma^2) \; , \; Cov\left(\mathbf{y}|\mathbf{X}\right) = \sigma^2\mathbf{I} $$
@ -150,11 +141,11 @@ Ainsi, un estimateur de la covariance est :
$$\hat{Cov}[\hat{\boldsymbol\beta}] = s^2 \left(\mathbf{X}^T\mathbf{X}\right)^{-1}$$
Les éléments diagonaux sont les carrés des erreurs standards pour les coefficients $\hat{\beta_i}$.
# Indicateur $R^2$ de la qualité d'une régression linéaire
## Indicateur $R^2$ de la qualité d'une régression linéaire
L'indicateur $R^2$ est le carré de la corrélation estimée entre $\hat{y}_i$ et $y_i$. Nous pouvons également dire que c'est le pourcentage de la variance de $y_i$ expliquée par $\mathbf{x_i}$. Il y a une version ajustée de cet estimateur qui est sinon naturellement biaisé vers le haut.
# Mesure des colinéarités par le facteur d'inflation de la variance
## Mesure des colinéarités par le facteur d'inflation de la variance
Pour mesurer les colinéarités entre la variable $\mathbf{x_1}$ et les variables $\mathbf{x_2},\dots,\mathbf{x_p}$, nous pouvons apprendre un modèle de la forme :
$$\mathbf{x_1} = \beta_0 + \beta_2\mathbf{x_2} + \dots + \beta_p\mathbf{x_p} + \epsilon$$
@ -162,7 +153,7 @@ Noton $R_1^2$ le carré de la corrélation estimée entre $\hat{\mathbf{x}}_1$ e
$$VIF_1 = \frac{1}{1 - R_1^2}$$
En pratique, quand $VIF(\beta_i)>10$, la colinéarité est importante.
# Effet levier des observations sur les prédictions
## Effet levier des observations sur les prédictions
$$
\begin{aligned}

2
02_moindres_carres.R

@ -6,4 +6,4 @@ polyreg2 <- function(data, degre) {
A = outer(xs, 0:degre, "^")
gram = t(A) %*% A
solve(gram, as.vector(t(A) %*% data$Y))
}
}

26
02_moindres_carres.Rmd

@ -1,18 +1,11 @@
---
title: "02 Méthode des moindres carrés"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Méthode des moindres carrés
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
source("02_moindres_carres.R", local = knitr::knit_global())
```
# Espace de fonctions
## Espace de fonctions
Soit un espace vectoriel composé de fonctions. Une base de cet espace est un ensemble de fonctions ($f_1, f_2, \dots f_p$) tel que toute fonction de l'espace s'exprime comme combinaison linéaire des fonctions de base.
\[
@ -27,16 +20,16 @@ f_1(\mathbf{x^{(2)}}) & f_2(\mathbf{x^{(2)}}) & \dots & f_p(\mathbf{x^{(2)}}) \\
f_1(\mathbf{x^{(n)}}) & f_2(\mathbf{x^{(n)}}) & \dots & f_p(\mathbf{x^{(n)}})
\end{array} \right)
\left( \begin{array}{c}
\beta_1 \\ \beta_2 \\ \dots \\ \beta_p
\beta_1 \\ \beta_2 \\ \dots \\ \beta_p
\end{array} \right)
=
\left( \begin{array}{c}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(n)}
\end{array} \right)
\]
Nous notons ce système linéaire $\mathbf{X}\boldsymbol\beta = \mathbf{y}$.
# Expression matricielle
## Expression matricielle
Le système linéaire $\mathbf{X}\boldsymbol\beta = \mathbf{y}$ avec $\mathbf{X} \in \mathbb{R}^{n \times p}$ n'a pas de solution quand le nombre d'observations dépasse le nombre de fonctions de base (c'est-à-dire, $n>p$). Une approche possible est alors de chercher une approximation $\mathbf{X}\boldsymbol\beta \approx \mathbf{y}$ qui minimise la somme des carrés des erreurs : $\|\mathbf{X}\boldsymbol\beta-\mathbf{y}\|^2_2$.
@ -67,7 +60,7 @@ Ainsi, quand $n>p$, la solution approximée $\hat{\boldsymbol\beta}$, telle que
\mathbf{X}^T\mathbf{X} \boldsymbol\beta = \mathbf{X}^T\mathbf{y}
\]
# Méthode des moindres carrés appliquée à la régression polynomiale
## Méthode des moindres carrés appliquée à la régression polynomiale
Pour un polynôme de degré $p-1$, les fonctions de base mentionnées ci-dessus sont : $f_1(x)=1$, $f_2(x)=x$, $f_3(x)=x^2$,..., $f_p(x)=x^{p-1}$. Elles permettent de définir la matrice des données $\mathbf{X}$ et la matrice de Gram $\mathbf{X}^T\mathbf{X}$.
@ -83,4 +76,9 @@ plt(data,f)
pltpoly(coef)
```
Ce polynôme de degré trois modélise mieux la fonction génératrice inconnue que celui de degré quatre qui ne commettait aucune erreur sur les données observées.
Ce polynôme de degré trois modélise mieux la fonction génératrice inconnue que celui de degré quatre qui ne commettait aucune erreur sur les données observées.
## Annexe code source
```{r, code=readLines("02_moindres_carres.R"), eval=FALSE}
```

2
03_tikhonov.R

@ -13,4 +13,4 @@ ridge <- function(alpha, data, degre) {
coef <- coef / attr(A,"scaled:scale")
inter <- Ym - coef %*% attr(A,"scaled:center")
coef <- c(inter, coef)
}
}

28
03_tikhonov.Rmd

@ -1,18 +1,11 @@
---
title: "03 Régularisation de Tikhonov"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Régularisation de Tikhonov
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
source("03_tikhonov.R", local = knitr::knit_global())
```
# Problèmes linéaires mal posés
## Problèmes linéaires mal posés
Avec moins d'observations que de fonctions de base ($n<p$), le système $\mathbf{X}\boldsymbol\beta=\mathbf{y}$ ne possède pas de solution unique. Même quand $n \geq p$, le système linéaire peut posséder une solution approchée préférable à la solution optimale. C'est en particulier vrai quand plusieurs observations sont très proches à un coefficient multiplicatif près (on parle de colinéarités). Par exemple, soit le système linéaire suivant :
\[
@ -30,7 +23,7 @@ Avec moins d'observations que de fonctions de base ($n<p$), le système $\mathbf
\]
Sa solution est $\boldsymbol\beta^T = (1001,-1000)$. Cependant, la solution approchée $\boldsymbol\beta^T = (0.5,0.5)$ semble préférable. En effet, la solution optimale a peu de chance de bien s'adapter à de nouvelles observations (par exemple, l'observation $(1,2)$ serait projetée sur l'étiquette $-999$).
# Ajout de contraintes de régularité
## Ajout de contraintes de régularité
Ainsi, lorsqu'il faut choisir entre plusieurs solutions, il peut être efficace d'exprimer une préférence envers celles dont les coefficients (ou paramètres) ont de faibles valeurs. Cela consiste par exemple à minimiser $|\beta_1|+|\beta_2|+\dots$ (aussi noté $\|\boldsymbol\beta\|_1$, la "norme 1") ou encore $\beta_1^2+\beta_2^2+\dots$ (aussi noté $\|\boldsymbol\beta\|_2^2$, le carré de la "norme 2"). Dans ce dernier cas, il s'agit de résoudre un nouveau problème de minimisation :
@ -49,7 +42,7 @@ C'est encore un problème de minimisation quadratique en $\boldsymbol\beta$ dont
En pratique, il s'agit donc d'ajouter une petite valeur positive $\lambda$ aux éléments de la diagonale de la matrice de Gram. Cette approche porte plusieurs noms dont "régularisation de Tikhonov" ou "régression Ridge".
# Standardisation
## Standardisation
La régression Ridge contraint les coefficients de la régression linéaire en les contractant vers 0. L'importance de la contraction est proportionnelle au carré des valeurs des coefficients. Ainsi, il est important que les domaines des différentes variables partage une même échelle. Sinon, l'intensité de l'influence d'une variable sur la régularisation dépendrait de son échelle.
@ -68,9 +61,9 @@ Nous pouvons également centrer ou standardiser la cible. Supposons que nous cho
Pour passer du modèle inféré sur les variables standardisées et la cible normalisée aux valeurs des paramètres sur les domaines originaux des variables, il faut donc ajouter une ordonnée à l'origine ("intercept") $\beta_0 = \mathbf{\bar{y}} - \sum_{j=1}^{p} \hat{\beta_j} \frac{\bar{x_j}}{\sigma_j}$ et remettre à l'échelle chaque paramètre $\beta_j$ en le multipliant par le facteur $\frac{1}{\sigma_j}$.
# Visulation de l'effet sur les paramètres de différents niveaux de régularisation
## Visulation de l'effet sur les paramètres de différents niveaux de régularisation
Sur notre exemple synthétique, nous affichons la fonction génératrice, le jeu de donnée et le polynôme de degré au plus sept découvert par régression ridge avec une valeur de $\lambda$ égale soit à $0$, soit à $10^{-4}$, soit à $1$.
Sur notre exemple synthétique, nous affichons la fonction génératrice, le jeu de données et le polynôme de degré au plus sept découvert par régression ridge avec une valeur de $\lambda$ égale soit à $0$, soit à $10^{-4}$, soit à $1$.
```{r}
set.seed(1123)
@ -104,7 +97,7 @@ matplot(lambdas, t(lcoef), type=c("b"), pch=1, col=1:8, log="x")
legend("topright", legend = 0:7, col=1:8, pch=1)
```
# Régularisation et complexité
## Régularisation et complexité
Essayons de mieux comprendre les raisons pour lesquelles la régularisation peut être efficace. Nous allons montrer qu'elle réduit la complexité d'un modèle prédictif. Dans ce contexte, qu'est-ce que la complexité ? C'est la sensibilité du modèle au bruit présent dans les données. Qu'est-ce que le bruit ? Il est possible de le définir après avoir fait l'hypothèse d'une famille de modèles qui auraient pu générer les données observées.
@ -122,4 +115,9 @@ Prenons l'exemple d'un modèle linéaire : $F(\mathbf{X}) = \mathbf{W}^T\mathbf{
& \|\mathbf{W}\|_2 \|\mathbf{\epsilon}\|_2 \\
\end{align*}
Donc, en pénalisant les grandes valeurs de $\|\mathbf{W}\|_2$, on réduit la sensibilité du modèle à de petites perturbations dans les données observées. Autrement dit, par l'ajout de régularisation, les prédictions associées aux plus proches voisins de $\mathbf{X}$ tendent à être similaires à celle associée à $\mathbf{X}$.
Donc, en pénalisant les grandes valeurs de $\|\mathbf{W}\|_2$, on réduit la sensibilité du modèle à de petites perturbations dans les données observées. Autrement dit, par l'ajout de régularisation, les prédictions associées aux plus proches voisins de $\mathbf{X}$ tendent à être similaires à celle associée à $\mathbf{X}$.
## Annexe code source
```{r, code=readLines("03_tikhonov.R"), eval=FALSE}
```

2
04_validation_croisee.R

@ -12,7 +12,7 @@ splitdata <- function(data,p) {
# lambdas[l] est une liste de valeurs pour l'hyperparamètre lambda.
# Notons Ridge[l] un modèle avec lambda <- lambdas[l].
# Découper aléatoirement le jeu de données d'entraînement en K plis F[i] disjoints.
# Découper aléatoirement le jeu d'entraînement en K plis F[i] disjoints.
# Pour l <- [1,...,L]
# Pour i <- [1,...,K]
# Apprendre Ridge[l] sur l'union des plis F[j] avec j!=i

18
04_validation_croisee.Rmd

@ -1,18 +1,11 @@
---
title: "04 Validation croisée"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Validation croisée
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
source("04_validation_croisee.R", local = knitr::knit_global())
```
# Principe de la validation croisée
## Principe de la validation croisée
Comment choisir la valeur du coefficient de régularisation $\lambda$ pour une régression ridge ? Notons en passant que $\lambda$ est un exemple de ce que l'on appelle un hyperparamètre car sa valeur doit être fixée avant de pouvoir apprendre les paramètres du modèle (dans notre cas, les coefficients d'un modèle linéaire).
@ -22,7 +15,7 @@ Il s'agirait donc de tester plusieurs valeurs de l'hyperparamètre $\lambda$ sur
Pour une approche souvent plus robuste, nous pouvons employer la stratégie dite de validation croisée à $K$ plis ("K Fold Cross-Validation"). Il s'agit de diviser aléatoirement le jeu d'entraînement en $K$ parties disjointes, appelées plis. Pour chaque valeur de l'hyperparamètre $\lambda$, nous apprenons $K$ modèles. Notons par exemple $M[\lambda_i,j]$, le j-ème des $K$ modèles appris pour la valeur $\lambda_i$ de l'hyperparamètre $\lambda$. Le jeu d'entraînement du modèle $M[\lambda_i,j]$ est constitué de l'union de $K-1$ plis, les $K$ plis initiaux auxquels on retire le j-ème pli qui joue le rôle de jeu de données de validation. Une estimation de l'erreur d'un modèle avec pour hyperparamètre $\lambda_i$ est obtenue en faisant la moyenne des erreurs des modèles $M[\lambda_i,j]$ sur les plis de validation. Ainsi, nous découvrons un meilleur hyperparamètre $\lambda_{best}$. Nous entraînons à nouveau un modèle sur l'ensemble du jeu d'entraînement (i.e., l'union des $K$ plis) avec un hyperparamètre $\lambda$ de valeur $\lambda_{best}$. Finalement, nous testons ce dernier modèle sur le jeu de test pour estimer son erreur sur des données encore jamais utilisées.
# Application de la validation croisée à la régularisation de Tikhonov
## Application de la validation croisée à la régularisation de Tikhonov
Le code source accompagnant ce chapitre comprend une fonction `splitdata` pour diviser le jeu de données en jeu d'entraînement et jeu de test. Ensuite, la fonction `kfoldridge` applique la stratégie de la validation croisée à $K$ plis sur le jeu d'entraînement pour une liste de valeurs de l'hyperparamètre $\lambda$. Elle retourne les coefficients du meilleur modèle et les moyennes des valeurs absolues des erreurs commises sur les plis de validation.
@ -59,3 +52,8 @@ Ce meilleur modèle atteint une erreur absolue moyenne de `r testmae` sur le jeu
plt(test,f)
pltpoly(reskfold$coef)
```
## Annexe code source
```{r, code=readLines("04_validation_croisee.R"), eval=FALSE}
```

156
05_b_svd_pca.Rmd

@ -1,20 +1,21 @@
---
title: "05-b SVD et Analyse en Composantes Principales"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# SVD et Analyse en Composantes Principales
# Projection sur les axes principaux
```{r, include=FALSE}
source("05_svd_pca.R", local = knitr::knit_global())
```
## SVD
```{r}
set.seed(1123)
```
## 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
### 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}$$
@ -23,25 +24,45 @@ 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
### 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
### 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,g}^2 = \mathcal{I} = \sum_\alpha \lambda_\alpha$.
## Contribution des variables aux axes principaux
### 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
## Implémentation
Nous écrivons une fonction `fa` (_factor analysis_) qui :
- utilise l'algorithme `k-means` pour découvrir une centaine de clusters à partir des données standardisées,
- applique une décomposition en valeurs singulières sur les centres standardisés des clusters,
- calcule :
- le pourcentage de variance expliquée par chaque axe principal (`prctPrcp`),
- les facteurs (`fact`), c'est-à-dire les coordonnées des observations sur les axes principaux,
- les contributions des observations aux axes principaux (`ctr`),
- les contributions des axes principaux aux écarts au centre d'inertie des observations (`cos2`),
- les contributions des variables aux axes principaux (`varctr`)
Nous écrivons une fonction `print.fa` qui affiche sur les axes principaux `d1` (par défaut $1$) et `d2` (par défaut $2$) les centres des clusters qui contribuent le plus à ces axes.
Nous écrivons une fonction `away.fa` qui retourne le cluster (son identifiant, sa taille et les noms des observations qui le composent) qui a le plus d'inertie (i.e., qui est le plus éloigné du centre d'inertie, c'est-à-dire l'origine du repère pour des données centrées) le long de l'axe principal `d` (par défaut $1$).
```{r, code=readLines("05_svd_pca.R"), eval=FALSE}
```
## 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}
@ -56,124 +77,77 @@ 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")]
X <- abalone[,c("length", "diameter", "height", "whole.wt","shucked.wt",
"viscera.wt", "shell.wt")]
```
Nous écrivons une fonction pour standardiser $\mathbf{X}$.
Nous réalisons une première fois l'analyse en composantes principales.
```{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)
}
fam <- fa(X) # fam pour 'factor analysis model'
```
```{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.
Nous affichons le pourcentage de variance expliquée par chaque axe principal.
```{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), "/")
fam$prctPrcp
```
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)
print(fam)
```
```{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])
far <- away(fam,2)
```
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 :
Le cluster `r far$id` (`far$id`) qui explique le plus la variance du deuxième axe principal semble anormalement important. Il contribue à expliquer `r round(fam$ctr[far$id,2]*100, 2)` % (`round(fam$ctr[far$id,2]*100, 2)`) de la variance du second axe. Il est composé de seulement `r far$size` (`far$size`) élément :
```{r}
abalone[f2.max.name,]
abalone[far$names,]
```
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)
abalone <- abalone[!(row.names(abalone) %in% far$names),]
X <- abalone[,c("length", "diameter", "height", "whole.wt","shucked.wt",
"viscera.wt", "shell.wt")]
fam <- fa(X)
```
Nous affichons le pourcentage de variance expliquée par chaque axe principal.
```{r}
fam$prctPrcp
```
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)
print(fam)
```
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,]
far <- away(fam,2)
```
Le cluster `r far$id` est intrigant. Sa taille est : `r far$size`. C'est encore un singleton qui correspond à l'observation `r far$names` du jeu de données original. Nous découvrons qu'il s'agit sans doute d'une anomalie pour la variable `height` :
```{r}
abalone[far$names,]
```
```{r}
boxplot(abalone$height)
```
Profitons-en pour calculer les contributions des axes principaux aux écarts au centre d'inertie des observations.
Nous affichons les contributions des axes principaux à l'écart au centre d'inertie du cluster `r far$id`.
```{r}
cos2 <- sweep(fact2, 1, rowSums(fact2), "/")
round(cos2[intrig,]*100,2)
fam$cos2[far$id,]
```
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)
fam$varctr
```
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.

2
05_b_svd_pca_slides.Rmd

@ -88,7 +88,7 @@ round((svd.C$d^2 / sum(svd.C$d^2))*100, 2)
# Exemple `abalone`
## Facteurs et des contributions des observations aux axes
## Facteurs et contributions des observations aux axes
- $\mathbf{F} \triangleq \mathbf{X}\mathbf{V} = \mathbf{U} \mathbf{D} \mathbf{V}^T\mathbf{V} = \mathbf{U} \mathbf{D}$
- $CTR_{i,\alpha} = \frac{f_{i,\alpha}^2}{\sum_i f_{i,\alpha}^2} = \frac{f_{i,\alpha}^2}{\lambda_\alpha}$

6
05_presentation_svd.R

@ -1,11 +1,11 @@
# chargement du package pixmap qui peut être installé avec la commande
# install.packages("pixmap");
library(pixmap) ;
library(pixmap) ;
read_grey_img <- function(file) {
img = read.pnm("images/hortensia.pgm");
# chaque entrée de img@grey est une intensité de gris comprise entre 0 et 1
img@grey;
img@grey;
}
compress_SVD <- function(X, nd) {
@ -17,4 +17,4 @@ print_grey_img <- function(Y,...) {
Y[Y<0]<-0; # après reconstruction approchée à partir du résultat du svd le
Y[Y>1]<-1; # domaine [0,1] de la matrice initiale peut ne plus être respecté
image(t(apply(Y,2,rev)), col=grey(seq(0,1,length=256)), axes=FALSE, ...)
}
}

22
05_presentation_svd.Rmd

@ -1,23 +1,16 @@
---
title: "05 Présentation de la décomposition en valeurs singulières (SVD)"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Présentation de la décomposition en valeurs singulières (SVD)
```{r, include=FALSE}
source("05_presentation_svd.R", local = knitr::knit_global())
```
# Contexte
## Contexte
Nous présentons la décomposition en valeurs singulières, une méthode générique de compression d'un tableau rectangulaire de données. Elle possède de nombreuses applications. En particulier, elle permet de calculer en un seul entraînement les régressions ridge pour toutes les valeurs possibles de l'hyperparamètre $\alpha$. Plus généralement, elle permet souvent d'adopter une perspective géométrique pour mieux comprendre les problèmes et algorithmes de l'apprentissage automatique à partir de données.
# Recherche d'un sous-espace qui minimise les carrés des écarts à l'espace d'origine
## Recherche d'un sous-espace qui minimise les carrés des écarts à l'espace d'origine
Soit un tableau de données $\mathbf{X}$ composé de $N$ observations en lignes, chacune décrite par $P$ variables en colonnes. $x_{ij}$ est la valeur de la variable $j$ pour l'observation $i$. Les vecteurs lignes forment $N$ points de l'espace $\mathbb{R}^P$. Les vecteurs colonnes forment $P$ points de l'espace $\mathbb{R}^N$. Nous cherchons à projeter le nuage des points lignes sur un sous-espace $\mathcal{H} \subset \mathbb{R}^P$, tout en minimisant les déformations.
Soit un tableau de données $\mathbf{X}$ composé de $N$ observations en lignes, chacune décrite par $P$ variables en colonnes. $x_{ij}$ est la valeur de la variable $j$ pour l'observation $i$. Les vecteurs lignes forment $N$ points de l'espace $\mathbb{R}^P$. Les vecteurs colonnes forment $P$ points de l'espace $\mathbb{R}^N$. Nous cherchons à projeter le nuage des points lignes sur un sous-espace $\mathcal{H} \subset \mathbb{R}^P$, tout en minimisant les déformations.
Pour commencer, nous considérons le meilleur sous-espace $\mathcal{H}$ à une dimension, c'est-à-dire une droite définie par son vecteur directeur unitaire $\mathbf{v}$ ("unitaire" signifie ici que sa norme euclidienne est égale à $1$, soit $\sqrt{\mathbf{v}^T\mathbf{v}}=1$, ou encore, $\mathbf{v}^T\mathbf{v}=1$). Soit $M_i$ un des $N$ points de $\mathbb{R}^P$. A ce point correspond le vecteur $\mathbf{OM_i}$ aussi noté $\mathbf{x_i}$ car ses coordonnées se lisent sur la i-ème ligne de $\mathbf{X}$. Soit $H_i$ la projection de $M_i$ sur la droite $\mathcal{H}$.
@ -122,7 +115,7 @@ Par ailleurs, nous obtenons une reconstitution approchée de rang $k$ de $\mathb
Considérons l'image d'un hortensia en fleurs.
![hortensia](images/hortensia.jpg)\
![hortensia](images/hortensia.jpg)
Nous convertissons cette image, originellement au format JPEG, en un format en niveaux de gris, simple à manipuler, appelé PGM (*Portable Grey Map*). Dans ce format non compressé, l'image est représentée par une matrice dont chaque valeur, comprise entre 0 et 1, représente le niveau de gris d'un pixel. Nous opérons cette transformation grâce au programme *imagemagick* avec la commande `convert hortensia.jpg hortensia.pgm`.
@ -139,3 +132,8 @@ print_grey_img(compress_SVD(X,64), main="d=64", asp=1);
print_grey_img(compress_SVD(X,128), main="d=128", asp=1);
print_grey_img(compress_SVD(X,256), main="d=256", asp=1);
```
## Annexe code source
```{r, code=readLines("05_presentation_svd.R"), eval=FALSE}
```

11
06_matrice_definie_non_negative.Rmd

@ -1,11 +1,4 @@
---
title: "06 Matrice définie non négative"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Matrice définie non négative
Avant de présenter le cœur de la preuve de l'existence de la décomposition en valeurs singulières, nous rappelons quelques propriétés des matrices symétriques définies non négatives. Une matrice définie non négative (souvent noté PSD pour "positive semi-definite") $\mathbf{M}$ représente une métrique car elle définit un produit scalaire et donc une géométrie, c'est-à-dire qu'elle donne un sens aux notions de distance et d'angle.
@ -53,7 +46,7 @@ Nous remarquons par exemple que l'équation d'un cercle selon la métrique $\mat
= \{ &\mathbf{M} = \mathbf{T}^T\mathbf{T}\} \\
& \mathbf{x}^T \mathbf{T}^T \mathbf{T} \mathbf{x} = c \\
= \{ &\text{Par définition du produit scalaire, on retrouve le cercle déformé en une ellipse par } \mathbf{T} \} \\
& \|\mathbf{T}\mathbf{x}\|_{\mathbf{I}}^2 = c \\
& \|\mathbf{T}\mathbf{x}\|_{\mathbf{I}}^2 = c \\
\end{align*}
Ainsi, une matrice symétrique définie non négative projette le cercle unité sur une ellipse dont les axes orthonormaux sont les vecteurs propres et les longueurs des axes sont les valeurs propres.

13
07_multiplicateurs_de_lagrange.Rmd

@ -1,11 +1,4 @@
---
title: "07 Optimisation sous contraintes et multiplicateurs de Lagrange"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Optimisation sous contraintes et multiplicateurs de Lagrange
Avant de pouvoir présenter la preuve de l'existence de la décomposition en valeurs singulières, nous devons introduire une stratégie souvent très utile pour résoudre un problème d'optimisation sous contraintes : la méthode dite des multiplicateurs de Lagrange.
@ -16,7 +9,7 @@ Sur cet exemple, les lignes de niveaux, ou contours, de $F$ sont des cercles cen
La solution $\mathbf{x*}$ appartient au contour de $F$ tangent à la contrainte $K$. Il faut donc que leurs gradients soient alignés : $\nabla F = \lambda \nabla K$.
\[
\nabla F =
\nabla F =
\left( \begin{array}{c}
\partial F / \partial x_1 \\
\partial F / \partial x_2
@ -27,7 +20,7 @@ La solution $\mathbf{x*}$ appartient au contour de $F$ tangent à la contrainte
2 x_2
\end{array} \right)
\quad ; \quad
\nabla K =
\nabla K =
\left( \begin{array}{c}
\partial K / \partial x_1 \\
\partial K / \partial x_2

11
08_derivation_svd.Rmd

@ -1,11 +1,4 @@
---
title: "08 Dérivation du SVD"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Dérivation du SVD
Après ces préambules sur la notion de métrique associée à une matrice définie non négative et sur la méthode d'optimisation dite des multiplicateurs de Lagrange, nous reprenons la dérivation de la preuve de l'existence de la décomposition en valeurs singulières pour toute matrice.
@ -13,7 +6,7 @@ Rappelons que nous avons noté $\mathbf{v_1}$ le vecteur directeur du meilleur s
Le résultat est sans difficulté plus général. Nous le prouvons pour toute matrice symétrique $\mathbf{A}$ (non seulement $\mathbf{X}^T\mathbf{X}$) et toute métrique de $\mathcal{R}^P$ représentée par une matrice symétrique définie non négative $\mathbf{M}$ (non seulement la matrice identité $\mathbf{I}$).
Nous rappelons que $\mathbf{v}^T\mathbf{A}\mathbf{v} = \Sigma a_{i,j}v_iv_j$. Ainsi, comme $\mathbf{A}$ et $\mathbf{M}$ sont symétriques (i.e., $a_{i,j}=a_{j,i}$ et $m_{i,j}=m_{j,i}$), les dérivées partielles des formes quadratiques ont les formes suivantes :
Nous rappelons que $\mathbf{v}^T\mathbf{A}\mathbf{v} = \Sigma a_{i,j}v_iv_j$. Ainsi, comme $\mathbf{A}$ et $\mathbf{M}$ sont symétriques (i.e., $a_{i,j}=a_{j,i}$ et $m_{i,j}=m_{j,i}$), les dérivées partielles des formes quadratiques ont les formes suivantes :
\[
\frac{\partial \mathbf{v}^T\mathbf{A}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{A}\mathbf{v} \quad \text{et} \quad \frac{\partial \mathbf{v}^T\mathbf{M}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{M}\mathbf{v}
\]

17
09_puissance_iteree_valeur_propre.Rmd

@ -1,15 +1,8 @@
---
title: "09 Plus grande valeur propre et puissance itérée"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Plus grande valeur propre et puissance itérée
Dans l'intention de comprendre plus tard comment calculer en pratique une décomposition en valeurs singulières, nous commençons par présenter la méthode dite de la puissance itérée pour calculer la plus grande valeur propre et son vecteur propre associé pour une matrice carrée.
# Décomposition en valeurs propres d'une matrice carrée
## Décomposition en valeurs propres d'une matrice carrée
Dire que $\mathbf{u}$ est un vecteur propre de la matrice $\mathbf{A}$ associé à la valeur propre $\lambda$ signifie que $\mathbf{A}\mathbf{u}=\lambda\mathbf{u}$. Remarquons que tout multiple de $\mathbf{u}$ est également un vecteur propre associé à la même valeur propre $\lambda$. Considérons l'exemple de la matrice ci-dessous.
\[
@ -30,7 +23,7 @@ Notons $\mathbf{e_1}^T=(1 \; 0 \; 0)$, $\mathbf{e_2}^T=(0 \; 1 \; 0)$ et $\mathb
Nous observons en particulier que la valeur propre $-1$ est associée à deux vecteurs propres linéairement indépendants : $\mathbf{e_3}$ et $(\mathbf{e_1}-\mathbf{e_2})$.
# Plus grande valeur propre et méthode de la puissance itérée
## Plus grande valeur propre et méthode de la puissance itérée
Posons $\mathbf{A}\in\mathbb{R}^{N\times N}$ et $\mathbf{u^{(0)}}\in\mathbb{R}^N$ avec $\|\mathbf{u^{(0)}}\|=1$. Supposons que $\mathbf{u^{(0)}}$ puisse se décomposer en une somme de vecteurs propres de $\mathbf{A}$.
@ -66,7 +59,7 @@ Par induction, nous obtenons :
\begin{align*}
& \mathbf{u^{(k)}} \\
= \{& \gamma_k = \| \lambda_1^k\mathbf{u_1} + \sum_{j=2}^{r} \lambda_j^k\mathbf{u_j} \| \} \\
= \{& \gamma_k = \| \lambda_1^k\mathbf{u_1} + \sum_{j=2}^{r} \lambda_j^k\mathbf{u_j} \| \} \\
& \frac{1}{\gamma_k} \left(\lambda_1^k\mathbf{u_1} + \sum_{j=2}^{r} \lambda_j^k\mathbf{u_j}\right) \\
= \phantom{\{}& \\
& \frac{\lambda_1^k}{\gamma_k} \left(\mathbf{u_1} + \sum_{j=2}^{r} \left(\frac{\lambda_j}{\lambda_1}\right)^k \mathbf{u_j}\right) \\
@ -80,7 +73,7 @@ Par ailleurs, $\mathbf{u^{(k)}}$ étant à chaque itération normalisé, nous av
& \frac{|\lambda_1^k|}{\gamma_k} = \frac{1}{\|\mathbf{u_1} + \sum_{j=2}^{r} \left(\frac{\lambda_j}{\lambda_1}\right)^k \mathbf{u_j}\|}
\end{align*}
D'où :
D'où :
\[
\lim\limits_{k \to \infty} \frac{|\lambda_1^k|}{\gamma_k} = \frac{1}{\|\mathbf{u_1}\|}
\]

11
10_puissance_iteree_svd.Rmd

@ -1,11 +1,4 @@
---
title: "10 Méthode de la puissance itérée et SVD"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Méthode de la puissance itérée et SVD
Nous montrons comment adapter la méthode de la puissance itérée au calcul de la plus grande valeur singulière et de ses deux vecteurs singuliers associés. Nous considérons une matrice de données rectangulaire $\mathbf{A}\in\mathbb{R}^{N\times M}$. Nous notons $\mathbf{U}\in\mathbb{R}^{N\times N}$ et $\mathbf{V}\in\mathbb{R}^{M\times M}$ les deux matrices orthogonales qui contiennent les vecteurs singuliers en colonnes. Nous notons enfin $\mathbf{\Sigma} = diag(\sigma_1,\sigma_2,\dots,\sigma_r,0,\dots)\in\mathbb{R}^{N\times M}$ la matrice qui contient sur sa diagonale principale les valeurs singulières, avec $r$ le range de $\mathbf{A}$. Nous avons donc :
@ -68,7 +61,7 @@ De même :
Donc :
\[
\frac{\|\mathbf{u^{(k+1)}}\|^2}{\|\mathbf{v^{(k)}}\|^2} = 1 =
\sigma_1^2 \left(\frac{\delta_{2k+1}^2}{\delta_{2k}^2}\right)
\sigma_1^2 \left(\frac{\delta_{2k+1}^2}{\delta_{2k}^2}\right)
\left(\frac{\sum_{j=1}^{n_1} y_j^2 \; + \; \sum_{j=n_1}^{r}\left(\frac{\sigma_j}{\sigma_1}\right)^{4k+2}y_j^2}
{\sum_{j=1}^{n_1} y_j^2 \; + \; \sum_{j=n_1}^{r}\left(\frac{\sigma_j}{\sigma_1}\right)^{4k}y_j^2}\right)
\]

17
11_projecteurs.Rmd

@ -1,15 +1,8 @@
---
title: "11 Projecteurs orthogonaux et SVD"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Projecteurs orthogonaux et SVD
Nous introduisons la notion de projecteur et nous montrons son lien avec la décomposition en valeurs signulières (SVD). Ainsi, nous comprendrons mieux la géométrie du SVD et nous serons plus tard en mesure de construire un algorithme qui est une extension de la méthode de la puissance itérée pour le calcul des $k$ premières valeurs singulières et vecteurs singuliers.
# Projecteur
## Projecteur
Un projecteur est une matrice carrée idempotente : $\mathbf{P}^2=\mathbf{P}$. Rappelons que l'espace des colonnes de $\mathbf{P}$ formé de l'ensemble des vecteurs qui s'expriment comme une combinaison linéaire des colonnes de $\mathbf{P}$, c'est-à-dire tous les vecteurs $\mathbf{v}$ qui s'écrivent sous la forme $\mathbf{P}\mathbf{x}$. L'espace des colonnes de $\mathbf{P}$ est également appelé l'image de $\mathbf{P}$ et se note $Im(\mathbf{P})$. Un vecteur $\mathbf{v}$ qui appartient à l'image d'un projecteur $\mathbf{P}$ est tel que $\mathbf{P}\mathbf{v}=\mathbf{P}$ :
@ -40,7 +33,7 @@ Le projecteur complémentaire à $\mathbf{P}$ est $\mathbf{I-P}$, il projette su
& \mathbf{v} \in Ker(\mathbf{P}) \\
= \{& \text{Par définition du noyau} \} \\
& \mathbf{P}\mathbf{v} = \mathbf{0} \\
\Rightarrow \phantom{\{}& \\
\Rightarrow \phantom{\{}& \\
& (\mathbf{I-P})\mathbf{v} = \mathbf{v} \\
= \{& \text{Par définition de l'image} \} \\
& \mathbf{v} \in Im(\mathbf{I-P}) \\
@ -50,7 +43,7 @@ Donc, $Ker(\mathbf{P}) \subseteq Im(\mathbf{I-P})$. Par ailleurs, un élément d
Donc, un projecteur $\mathbf{P}$ sépare l'espace initial en deux sous-espaces $S1 = Im(\mathbf{P})$ et $S2 = Ker(\mathbf{P})$. $\mathbf{P}$ projette **sur** $S1$ **parallèlement à** (ou **le long de**) $S2$.
# Projecteur orthogonal
## Projecteur orthogonal
Un projecteur est dit orthogonal quand $S1$ et $S2$ sont orthogonaux. Montrons qu'une projection $\mathbf{P}$ est orthogonale si et seulement si $\mathbf{P}=\mathbf{P}^T$. Montrons d'abord que $\mathbf{P}=\mathbf{P}^T \Rightarrow \text{orthogonalité}$.
@ -110,7 +103,7 @@ Il suffit donc de montrer que chaque $\mathbf{A}\mathbf{v_i}$ est nul.
& 0 \\
\end{align*}
# Projection sur une base quelconque
## Projection sur une base quelconque
Nous cherchons à projeter un vecteur $\mathbf{v}$ sur l'espace des colonnes d'une matrice $\mathbf{A}$. Notons le résultat de cette projection $\mathbf{y} \in Im(\mathbf{A})$. $\mathbf{y}-\mathbf{v}$ est orthogonal à $Im(\mathbf{A})$ (c'est-à-dire qu'il appartient à $Ker(\mathbf{A})$). Autrement dit, en posant $\mathbf{y}=\mathbf{A}\mathbf{x}$, nous avons :

64
12_factorisation_qr.Rmd

@ -1,19 +1,10 @@
---
title: "12 Factorisation QR"
output:
bookdown::pdf_document2:
number_section: yes
extra_dependencies:
algorithm2e: [ruled,vlined,linesnumbered]
toc: false
classoption: fleqn
---
# Factorisation QR
```{r, include=FALSE}
source("12_factorisation_qr.R", local = knitr::knit_global())
```
# Factorisation QR
## Factorisation QR
Nous introduisons une décomposition matricielle qui nous permettra plus tard d'implémenter les calculs des décompositions en valeurs propres et en valeurs singulières.
@ -27,7 +18,7 @@ Nous cherchons des vecteurs orthonormaux qui forment une base pour les espaces c
Nous avons $<\mathbf{a_1}> \subseteq <\mathbf{a_1}, \mathbf{a_2}> \subseteq \dots \subseteq <\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_N}>$.
Nous cherchons des vecteurs orthogonaux $\mathbf{q_i}$ tels que :
Nous cherchons des vecteurs orthogonaux $\mathbf{q_i}$ tels que :
\[
<\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_j}> = <\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_j}> \quad \text{pour } j=1,2,\dots,N
\]
@ -48,7 +39,7 @@ Sans perte de généralité, considérons la situation où $M>N$. Nous appelons
Nous pouvons aussi considérer la forme dite complète de la factorisation QR, avec $\mathbf{Q} \in \mathbb{R}^{M \times M}$ et $\mathbf{R} \in \mathbb{R}^{M \times N}$. Dans ce cas, les $M-N$ dernières lignes de $\mathbf{R}$ sont nulles et il suffit de compléter $\mathbf{\hat{Q}}$ en lui ajoutant $M-N$ colonnes orthogonales.
# Orthogonalisation de Gram-Schmidt
## Orthogonalisation de Gram-Schmidt
Nous présentons l'algorithme de Gram-Schmidt pour la décomposition QR.
@ -92,24 +83,24 @@ telles que $\mathbf{A} = \mathbf{Q} \mathbf{R}$}
$\mathbf{q_j} \gets \mathbf{v_j} / r_{jj}$ \;
}
\caption{Gram-Schmidt classique (GS)}
\end{algorithm}
\end{algorithm}
# Algorithme de Gram-Schmidt modifié
## Algorithme de Gram-Schmidt modifié
## Dérivation de l'algorithme de Gram-Schmidt modifié
### Dérivation de l'algorithme de Gram-Schmidt modifié
Nous présentons une autre manière de voir ce résultat qui nous mènera à un algorithme équivalent mais numériquement plus stable.
Nous pouvons considérer que $\mathbf{q_j}$ est obtenu en projetant $\mathbf{a_j}$ sur l'espace orthogonal à $<\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_{j-1}}>$, cet espace étant par construction égal à $<\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}>$ :
\[
\mathbf{q_j} = \frac{\mathbf{P_j}\mathbf{a_j}}{\|\mathbf{P_j}\mathbf{a_j}\|} \quad ; \quad
\mathbf{P_j} = \mathbf{I} - \mathbf{Q_{j-1}}\mathbf{Q_{j-1}}^T \; \text{avec }
\mathbf{P_j} = \mathbf{I} - \mathbf{Q_{j-1}}\mathbf{Q_{j-1}}^T \; \text{avec }
\mathbf{Q_{j-1}} = (\mathbf{q_1} \; \mathbf{q_2} \; \dots \mathbf{q_{j-1}})
\]
De façon équivalente, $\mathbf{P_j}$ peut s'exprimer comme le produit de projecteurs sur les espaces orthogonaux à $\mathbf{q_1}$, puis $\mathbf{q_2}$, etc. jusqu'à $\mathbf{q_{j-1}}$ :
\[
\mathbf{P_j} = \mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots
\mathbf{P_j} = \mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots
\mathbf{P_{\bot \mathbf{q_{2}}}} \mathbf{P_{\bot \mathbf{q_{1}}}} \quad
\text{avec } \mathbf{P_{\bot \mathbf{q}}} = \mathbf{I} - \mathbf{q}\mathbf{q}^T
\]
@ -118,7 +109,7 @@ Dans la version classique de Gram-Schmidt, $\mathbf{v_j}=\mathbf{P_j}\mathbf{a_j
Dans la version modifiée de Gram-Schmidt :
\[
\mathbf{v_j}=\mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots
\mathbf{v_j}=\mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots
\mathbf{P_{\bot \mathbf{q_{2}}}} \mathbf{P_{\bot \mathbf{q_{1}}}}
\mathbf{a_j}
\]
@ -153,7 +144,7 @@ telles que $\mathbf{A} = \mathbf{Q} \mathbf{R}$}
\caption{Gram-Schmidt modifié (MGS)}
\end{algorithm}
## Test de l'algorithme
### Test de l'algorithme
Nous comparons notre implémentation (voir la fonction `mgs` dans le code source associé à ce module) avec le résultat d'un petit exemple issu de [Wikipedia](https://en.wikipedia.org/wiki/Gram–Schmidt_process#Example).
@ -163,7 +154,7 @@ qr <- mgs(mat)
all.equal(qr$Q, (1/(sqrt(10))) * matrix(c(3,1,-1,3),nrow=2,ncol=2))
```
## Comparaison de la stabilité numérique des deux versions de l'algorithme de Gram-Schmidt
### Comparaison de la stabilité numérique des deux versions de l'algorithme de Gram-Schmidt
Observons comment de petites erreurs de calculs se propagent dans le cas de l'algorithme de Gram-Schmidt (GS). Mettons qu'une petite erreur ait été commise sur le calcul de $\mathbf{q_2}$ qui s'en trouve de ce fait n'être pas parfaitement perpendiculaire à $\mathbf{q_1}$ :
\[
@ -173,11 +164,11 @@ Observons comment de petites erreurs de calculs se propagent dans le cas de l'al
Quelle est l'effet de cette erreur sur la valeur de $\mathbf{v_3}$ ?
\begin{align*}
\mathbf{v_3} &= \mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\mathbf{q_1} -
\mathbf{v_3} &= \mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\mathbf{q_1} -
(\mathbf{q_2}^T\mathbf{a_3})\mathbf{q_2} \\
\mathbf{q_2}^T\mathbf{v_3} &= \mathbf{q_2}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\epsilon -
\mathbf{q_2}^T\mathbf{v_3} &= \mathbf{q_2}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\epsilon -
(\mathbf{q_2}^T\mathbf{a_3}) = - (\mathbf{q_1}^T\mathbf{a_3})\epsilon \\
\mathbf{q_1}^T\mathbf{v_3} &= \mathbf{q_1}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3}) -
\mathbf{q_1}^T\mathbf{v_3} &= \mathbf{q_1}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3}) -
(\mathbf{q_2}^T\mathbf{a_3})\epsilon = - (\mathbf{q_2}^T\mathbf{a_3})\epsilon \\
\end{align*}
@ -198,14 +189,33 @@ Il n'y a pas d'erreur pour l'orthogonalité entre $\mathbf{q_3}$ et $\mathbf{q_2
\mathbf{q_1}^T\mathbf{v_3} &= \mathbf{q_1}^T\mathbf{v_3}^{(1)} - \mathbf{q_2}^T\mathbf{v_3}^{(1)}\epsilon \\
\mathbf{q_2}^T\mathbf{v_3}^{(1)} &= \mathbf{q_2}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon \\
\mathbf{q_1}^T\mathbf{v_3}^{(1)} &= \mathbf{q_1}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)} = 0 \\
\mathbf{q_1}^T\mathbf{v_3} &= -(\mathbf{q_2}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon) =
\mathbf{q_1}^T\mathbf{v_3} &= -(\mathbf{q_2}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon) =
-\mathbf{q_2}^T\mathbf{v_3}^{(0)}\epsilon + \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon^2
\end{align*}
Ainsi, comme pour la version classique de Gram-Schmidt, une erreur de l'ordre de grandeur de $\epsilon$ est commise sur l'orthogonalité entre $\mathbf{q_1}$ et $\mathbf{q_3}$. Par contre, il n'y a pas de propagation de l'erreur pour l'orthogonalité entre $\mathbf{q_2}$ et $\mathbf{q_3}$. C'est pourquoi l'algorithme de Gram-Schmidt modifié est numériquement plus stable que l'algorithme de Gram-Schmidt classique.
# Autres algorithmes
## Autres algorithmes
Notons qu'une autre approche très souvent employée pour le calcul de la décomposition QR est la méthode dite de Householder. Elle assure le maintien d'une quasi parfaite orthogonalité des colonnes de $\mathbf{Q}$. Il existe également l'approche dite des rotations de Givens, plus facilement parallélisable et permettant la mise à jour de la décomposition suite à l'ajout d'une ligne à la matrice $\mathbf{X}$.
Par contre, Gram-Schmidt permet de calculer la forme réduite de la décomposition QR, alors que les algorithmes de Givens et Householder calculent la forme complète. Ainsi, Gram-Schmidt n'est pas tout à fait aussi stable mais peut parfois être plus efficace.
Par contre, Gram-Schmidt permet de calculer la forme réduite de la décomposition QR, alors que les algorithmes de Givens et Householder calculent la forme complète. Ainsi, Gram-Schmidt n'est pas tout à fait aussi stable mais peut parfois être plus efficace.
## Décomposition QR pour résoudre un système linéaire au sens des moindres carrés
Une transformation orthogonale $\mathbf{Q}$ ne change pas la norme du vecteur transformé.
$$
\|\mathbf{Q} \mathbf{y}\|_2^2 = \left(\mathbf{Q} \mathbf{y}\right)^T\left(\mathbf{Q} \mathbf{y}\right) = \mathbf{y}^T(\mathbf{Q}^T(\mathbf{Q}\mathbf{y} = \mathbf{y}^T\mathbf{y} = \|\mathbf{y}\|_2^2
$$
Avec $\mathbf{A}=\mathbf{Q}\mathbf{R}$ (décomposition QR en forme réduite, avec $\mathbf{Q} \in \mathbb{R}^{M \times N}$ et $\mathbf{R} \in \mathbb{R}^{N \times N}$) :
$$
\| \mathbf{A}\mathbf{x} - \mathbf{b}\|_2 = \| \mathbf{Q}^T \left(\mathbf{A}\mathbf{x} - \mathbf{b}\right)\|_2 = \| \mathbf{Q}^T \left(\mathbf{Q}\mathbf{R}\mathbf{x} - \mathbf{b}\right)\|_2 = \| \mathbf{R}\mathbf{x} - \mathbf{Q}^T\mathbf{b}\|_2
$$
Or, si la matrice $\mathbf{R}\mathbf{x}$ n'est pas singulière (i.e., s'il n'y a pas de colinéarités entre ses colonnes) le système linéaire $\mathbf{R}\mathbf{x} = \mathbf{Q}^T\mathbf{b}$, à $N$ équations et $N$ inconnues, se résout simplement par substitutions car $\mathbf{R}$ est de forme diagonale supérieure.
## Annexe code source
```{r, code=readLines("12_factorisation_qr.R"), eval=FALSE}
```

17
13_puissance_iteree_par_blocs.Rmd

@ -1,15 +1,6 @@
---
title: "13 Puissance itérée par bloc et SVD"
output:
bookdown::pdf_document2:
number_section: yes
extra_dependencies:
algorithm2e: [ruled,vlined,linesnumbered]
toc: false
classoption: fleqn
---
# Puissance itérée par bloc et SVD
# Principe et application à la décomposition en valeurs propres
## Principe et application à la décomposition en valeurs propres
Nous avons vu précédemment que la méthode de la puissance itérée permet de calculer la plus grande valeur propre et un vecteur propre associé pour une matrice carrée $\mathbf{A}$. La décomposition QR est au cœur d'une approche élégante pour étendre la méthode de la puissance itérée au calcul des $S$ plus grandes valeurs propres et leurs vecteurs propres associés.
@ -39,9 +30,9 @@ Les valeurs propres correspondantes sont sur la diagonale de $\mathbf{\Lambda}$.
$err \gets \|\mathbf{A}\mathbf{V} - \mathbf{V}\mathbf{\Lambda}\|$ \;
}
\caption{Puissance itérée par bloc pour les vecteurs propres}
\end{algorithm}
\end{algorithm}
# Application à la décomposition en valeurs singulières
## Application à la décomposition en valeurs singulières
Nous pouvons de même adapter l'algorithme de la puissance itérée au cas du SVD pour calculer les $S$ plus grandes valeurs singulières et les vecteurs singuliers associés pour une matrice $\mathbf{A} \in \mathbb{R}^{M \times N}$.

4
14_geometrie_ridge_svd.R

@ -50,7 +50,7 @@ ridge.svd <- function(data, degre, fold = FALSE) {
# alphas[l] est une liste de valeurs pour l'hyperparamètre alpha.
# Notons Ridge[l] un modèle avec alpha <- alphas[l].
# Découper aléatoirement le jeu de données d'entraînement en K plis F[i] disjoints.
# Découper aléatoirement le jeu d'entraînement en K plis F[i] disjoints.
# Pour l <- [1,...,L]
# Pour i <- [1,...,K]
# Apprendre Ridge[l] sur l'union des plis F[j] avec j!=i
@ -85,4 +85,4 @@ kfoldridge <- function(K, alphas, data, degre) {
ridge <- ridge.svd(data, degre)
coef <- ridge(bestalpha)
list(coef = coef, maes = maes, alpha = bestalpha)
}
}

34
14_geometrie_ridge_svd.Rmd

@ -1,20 +1,11 @@
---
title: "14 Géométrie de la régression ridge et SVD"
output:
bookdown::pdf_document2:
number_section: yes
extra_dependencies:
algorithm2e: [ruled,vlined,linesnumbered]
toc: false
classoption: fleqn
---
# Géométrie de la régression ridge et SVD
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
source("14_geometrie_ridge_svd.R", local = knitr::knit_global())
```
# Coefficients de la régression ridge en fonction du SVD
## Coefficients de la régression ridge en fonction du SVD
Exprimons le calcul des coefficients d'une régression ridge en fonction de la décomposition en valeurs singulières de la matrice des données observées $\mathbf{X} \in \mathbb{R}^{M \times N}$.
@ -22,23 +13,23 @@ Exprimons le calcul des coefficients d'une régression ridge en fonction de la d
& \mathbf{\hat{\beta}_\lambda} \\
= \{& \text{Voir la dérivation de la régression ridge dans un précédent module.} \} \\
& (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y} \\
= \{& \text{SVD de $\mathbf{X}$ : } \mathbf{U}\mathbf{D}\mathbf{V}^T \quad
\mathbf{U} \in \mathbb{R}^{M \times M} \; , \;
= \{& \text{SVD de $\mathbf{X}$ : } \mathbf{U}\mathbf{D}\mathbf{V}^T \quad
\mathbf{U} \in \mathbb{R}^{M \times M} \; , \;
\mathbf{D} \in \mathbb{R}^{M \times N} \; , \;
\mathbf{V} \in \mathbb{R}^{N \times N} \} \\
& (\mathbf{V}\mathbf{D}^T\mathbf{U}^T\mathbf{U}\mathbf{D}\mathbf{V}^T + \lambda\mathbf{I})^{-1}
& (\mathbf{V}\mathbf{D}^T\mathbf{U}^T\mathbf{U}\mathbf{D}\mathbf{V}^T + \lambda\mathbf{I})^{-1}
\mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\
= \{& \text{$\mathbf{U}$ et $\mathbf{V}$ sont orthogonales : } \mathbf{I} = \mathbf{V}\mathbf{V}^T \} \\
& (\mathbf{V}\mathbf{D}^T\mathbf{D}\mathbf{V}^T + \lambda\mathbf{V}\mathbf{V}^T)^{-1} \mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\
= \phantom{\{}& \\
& (\mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})\mathbf{V}^T)^{-1} \mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\
= \{& (\mathbf{X}\mathbf{Y})^{-1} = \mathbf{Y}^{-1} \mathbf{X}^{-1} \quad
= \{& (\mathbf{X}\mathbf{Y})^{-1} = \mathbf{Y}^{-1} \mathbf{X}^{-1} \quad
\text{$\mathbf{V}$ est orthogonale : } \mathbf{V}^{-1} = \mathbf{V}^T\} \\
& \mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})^{-1}\mathbf{V}^T\mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\
= \phantom{\{}& \\
& \mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})^{-1}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\
= \{& \text{Soit $d_j$ le jème élément sur la diagonale de $\mathbf{D}$,
$\mathbf{u_j}$ et $\mathbf{v_j}$ les jèmes colonnes de respectivement $\mathbf{U}$ et $\mathbf{V}$} \} \\
= \{& \text{Soit $d_j$ le jème élément sur la diagonale de $\mathbf{D}$,
$\mathbf{u_j}$ et $\mathbf{v_j}$ les jèmes colonnes de resp. $\mathbf{U}$ et $\mathbf{V}$} \} \\
& \sum_{d_j>0} \mathbf{v_j} \frac{d_j}{d_j^2 + \lambda} \mathbf{u_j}^T\mathbf{y} \\
\end{align*}
@ -102,7 +93,7 @@ plt(test,f)
pltpoly(reskfold$coef)
```
# Régression ridge et géométrie
## Régression ridge et géométrie
Observons la relation entre les étiquettes prédites $\mathbf{\hat{y}_\lambda}$ et les étiquettes observées $\mathbf{y}$.
@ -114,4 +105,9 @@ Observons la relation entre les étiquettes prédites $\mathbf{\hat{y}_\lambda}$
Nous remarquons qu'en labsence de régularisation, $\lambda = 0$, les valeurs estimées $\mathbf{\hat{y}}$ sont les projections sur les axes principaux $\mathbf{u_j}$ -- qui couvrent l'espace des colonnes de $\mathbf{X}$, i.e. $Im(\mathbf{X})$ -- des valeurs observées $\mathbf{y}$.
En présence de régularisation, $\lambda > 0$, les coordonnées, sur les axes principaux, de l'estimation $\mathbf{\hat{y}_\lambda}$ sont de plus en plus contractées lorsqu'on progresse vers les axes qui expliquent de moins en moins la variabilités des données.
En présence de régularisation, $\lambda > 0$, les coordonnées, sur les axes principaux, de l'estimation $\mathbf{\hat{y}_\lambda}$ sont de plus en plus contractées lorsqu'on progresse vers les axes qui expliquent de moins en moins la variabilités des données.
## Annexe code source
```{r, code=readLines("14_geometrie_ridge_svd.R"), eval=FALSE}
```

62
15_loocv.R

@ -1,29 +1,37 @@
# 15 LOOCV
ridge.svd <- function(data, degre) {
X <- data$X
ridgeSVD <-
function(X, y)
{
n <- nrow(X)
A <- scale(outer(c(X), 1:degre, "^"))
Ym <- mean(data$Y)
Y <- data$Y - Ym
As <- svd(A)
d <- As$d
X.init <- X
y.init <- y
X <- scale(X)
ym <- mean(y)
y <- y - ym
Xs <- svd(X)
d <- Xs$d
function(lambda) {
coef <- c(As$v %*% ((d / (d^2 + lambda)) * (t(As$u) %*% Y)))
coef <- coef / attr(A,"scaled:scale")
inter <- Ym - coef %*% attr(A,"scaled:center")
coef <- c(Xs$v %*% ((d / (d^2 + lambda)) * (t(Xs$u) %*% y)))
coef <- coef / attr(X,"scaled:scale")
inter <- ym - coef %*% attr(X,"scaled:center")
coef <- c(inter, coef)
trace.H <- sum(d^2 / (d^2 + lambda))
pred <- polyeval(coef, X)
gcv <- sum( ((Y - pred) / (1 - (trace.H / n))) ^ 2 ) / n
yh <- coef[1] + X.init%*%coef[-1]
gcv <- sum( ((y.init - yh) / (1 - (trace.H / n))) ^ 2 ) / n
list(coef = coef, gcv = gcv)
}
}
ridge.loocv <- function(data, deg, lambdas) {
ridge <-
function(X, y, lambdas=NULL)
{
X <- as.matrix(X)
p <- ncol(X)
if(is.null(lambdas)) { lambdas <- 10^seq(-8,8,by=0.5) }
errs <- double(length(lambdas))
coefs <- matrix(data = NA, nrow = length(lambdas), ncol = 1+deg)
ridge <- ridge.svd(data, deg)
coefs <- matrix(data = NA, nrow = length(lambdas), ncol = p+1)
ridge <- ridgeSVD(X, y)
idx <- 1
for(lambda in lambdas) {
res <- ridge(lambda)
@ -34,7 +42,23 @@ ridge.loocv <- function(data, deg, lambdas) {
err.min <- min(errs)
lambda.best <- lambdas[which(errs == err.min)]
coef.best <- coefs[which(errs == err.min),]
pred <- polyeval(coef.best, data$X)
mae <- mean(abs(pred - data$Y))
list(coef = coef.best, lambda = lambda.best, mae = mae)
}
yh <- coef.best[1] + X%*%coef.best[-1]
mae <- mean(abs(yh - y))
r <- list(coef = coef.best,
lambda = lambda.best,
mae = mae,
coefs=coefs,
lambdas=lambdas)
class(r) <- "ridge"
return(r)
}
predict.ridge <-
function(o, newdata)
{
newdata <- as.matrix(newdata)
if(length(o$coef)-1!=ncol(newdata)) {
stop("Not the same number of variables btwn fitted ridge object and new data")
}
yh <- o$coef[1] + newdata%*%o$coef[-1]
}

27
15_loocv.Rmd

@ -1,11 +1,4 @@
---
title: "15 Validation croisée un contre tous"
output:
bookdown::pdf_document2:
number_section: yes
toc: false
classoption: fleqn
---
# Validation croisée un contre tous
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
@ -142,16 +135,17 @@ data = gendat(n,0.2)
splitres <- splitdata(data,0.8)
entr <- splitres$entr
test <- splitres$test
lambdas <- c(1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
ridge.res <- ridge.loocv(entr, deg, lambdas)
lambdas <- 10^seq(2,-8,by=-1)
entr.poly <- outer(c(entr$X), 1:deg, "^")
rm <- ridge(entr.poly, entr$Y, lambdas)
plt(entr,f)
pltpoly(ridge.res$coef)
pltpoly(rm$coef)
```
Ci-dessus, nous avons généré un jeu de données composé de `r n` observations et nous avons calculé par validation croisée un contre tous un polynôme de degré au plus égal à `r deg` qui modélise au mieux ces données. La valeur de $\lambda$ retenue est : `r ridge.res$lambda` et l'erreur absolue moyenne sur le jeu d'entraînement est : `r ridge.res$mae`.
Ci-dessus, nous avons généré un jeu de données composé de `r n` observations et nous avons calculé par validation croisée un contre tous un polynôme de degré au plus égal à `r deg` qui modélise au mieux ces données. La valeur de $\lambda$ retenue est : `r rm$lambda` et l'erreur absolue moyenne sur le jeu d'entraînement est : `r rm$mae`.
```{r}
testpred <- polyeval(ridge.res$coef, test$X)
testpred <- polyeval(rm$coef, test$X)
testmae <- mean(abs(testpred - test$Y))
```
@ -159,5 +153,10 @@ Ce meilleur modèle atteint une erreur absolue moyenne de `r testmae` sur le jeu
```{r}
plt(test,f)
pltpoly(ridge.res$coef)
pltpoly(rm$coef)
```
## Annexe code source
```{r, code=readLines("15_loocv.R"), eval=FALSE}
```

63
15_loocv_slides.Rmd

@ -57,22 +57,20 @@ $$
\small
```{r, eval=FALSE}
ridge.svd <- function(data, degre) {
X <- data$X
n <- nrow(X)
A <- scale(outer(c(X), 1:degre, "^"))
Ym <- mean(data$Y)
Y <- data$Y - Ym
As <- svd(A)
d <- As$d
ridgeSVD <-
function(X, y)
{
n <- nrow(X); X.init <- X; y.init <- y
X <- scale(X); ym <- mean(y); y <- y - ym
Xs <- svd(X); d <- Xs$d
function(lambda) {
coef <- c(As$v %*% ((d / (d^2 + lambda)) * (t(As$u) %*% Y)))
coef <- coef / attr(A,"scaled:scale")
inter <- Ym - coef %*% attr(A,"scaled:center")
coef <- c(Xs$v %*% ((d / (d^2 + lambda)) * (t(Xs$u) %*% y)))
coef <- coef / attr(X,"scaled:scale")
inter <- ym - coef %*% attr(X,"scaled:center")
coef <- c(inter, coef)
trace.H <- sum(d^2 / (d^2 + lambda))
pred <- polyeval(coef, X)
gcv <- sum( ((Y - pred) / (1 - (trace.H / n))) ^ 2 ) / n
yh <- coef[1] + X.init%*%coef[-1]
gcv <- sum( ((y.init - yh) / (1 - (trace.H / n))) ^ 2 ) / n
list(coef = coef, gcv = gcv)
}
}
@ -83,11 +81,10 @@ ridge.svd <- function(data, degre) {
\small
```{r, eval=FALSE}
ridge.loocv <- function(data, deg, lambdas) {
errs <- double(length(lambdas))
coefs <- matrix(data = NA, nrow = length(lambdas), ncol = 1+deg)
ridge <- ridge.svd(data, deg)
idx <- 1
ridge <- function(X, y, lambdas) {
p <- ncol(X); errs <- double(length(lambdas))
coefs <- matrix(data = NA, nrow = length(lambdas), ncol = p+1)
ridge <- ridgeSVD(X, y); idx <- 1
for(lambda in lambdas) {
res <- ridge(lambda)
coefs[idx,] <- res$coef
@ -97,9 +94,11 @@ ridge.loocv <- function(data, deg, lambdas) {
err.min <- min(errs)
lambda.best <- lambdas[which(errs == err.min)]
coef.best <- coefs[which(errs == err.min),]
pred <- polyeval(coef.best, data$X)
mae <- mean(abs(pred - data$Y))
list(coef = coef.best, lambda = lambda.best, mae = mae)
yh <- coef.best[1] + X%*%coef.best[-1]
mae <- mean(abs(yh - y))
r <- list(coef = coef.best, lambda = lambda.best, mae = mae,
coefs=coefs, lambdas=lambdas)
class(r) <- "ridge" ; return(r)
}
```
\normalsize
@ -115,10 +114,11 @@ data = gendat(n,0.2)
splitres <- splitdata(data,0.8)
entr <- splitres$entr
test <- splitres$test
lambdas <- c(1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
ridge.res <- ridge.loocv(entr, deg, lambdas)
lambdas <- 10^seq(2,-8,by=-1)
entr.poly <- outer(c(entr$X), 1:deg, "^")
rm <- ridge(entr.poly, entr$Y, lambdas)
plt(entr,f)
pltpoly(ridge.res$coef)
pltpoly(rm$coef)
```
\normalsize
@ -132,15 +132,16 @@ data = gendat(n,0.2)
splitres <- splitdata(data,0.8)
entr <- splitres$entr
test <- splitres$test
lambdas <- c(1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
ridge.res <- ridge.loocv(entr, deg, lambdas)
lambdas <- 10^seq(2,-8,by=-1)
entr.poly <- outer(c(entr$X), 1:deg, "^")
rm <- ridge(entr.poly, entr$Y, lambdas)
plt(entr,f)
pltpoly(ridge.res$coef)
pltpoly(rm$coef)
```
```{r, echo=FALSE}
bestLambda <- ridge.res$lambda
maeTrain <- ridge.res$mae
bestLambda <- rm$lambda
maeTrain <- rm$mae
```
- $\lambda$ : `r bestLambda`
@ -150,13 +151,13 @@ maeTrain <- ridge.res$mae
```{r, echo=FALSE}
plt(test,f)
pltpoly(ridge.res$coef)
pltpoly(rm$coef)
```
# Exemple - jeu de test
```{r}
testpred <- polyeval(ridge.res$coef, test$X)
testpred <- polyeval(rm$coef, test$X)
testmae <- mean(abs(testpred - test$Y))
```

84
16_biais_variance_estimateur.Rmd

@ -1,15 +1,6 @@
---
title: "16 Biais et variance d'un estimateur"
output:
bookdown::pdf_document2:
number_section: yes
includes:
in_header: preamble.tex
toc: false
classoption: fleqn
---
# Biais et variance d'un estimateur
# Biais et variance pour l'estimateur d'un paramètre d'un modèle paramétrique
## Biais et variance pour l'estimateur d'un paramètre d'un modèle paramétrique
Nous notons $\hat{\mu}$ la moyenne ($\hat{\mu} = E[\mathbf{x}]$) et $\hat{\sigma}^2$ la variance ($\hat{\sigma}^2 = E[(\mathbf{x} - E[\mathbf{x}])(\mathbf{x} - E[\mathbf{x}])]$) d'un modèle paramétrique $P(\mathbf{x};\mu,\sigma)$ dont on fait l'hypothèse qu'il aurait généré des données observées $\mathbf{x_1}, \mathbf{x_2},\dots, \mathbf{x_n}$. Étant données les observations, nous nous intéressons à des estimateurs $\hat{\mu}_{est}$ et $\hat{\sigma}_{est}$ des paramètres $\hat{\mu}$ et $\hat{\sigma}$.
@ -17,7 +8,7 @@ Un estimateur est non biaisé si sa moyenne sur l'ensemble de tous les jeux de d
En plus d'un faible biais, un bon estimateur possède une faible variance : $Var(\hat{\mu}_{est}) = E[(\hat{\mu} - \hat{\mu}_{est})^2]$, $Var(\hat{\sigma}_{est}) = E[(\hat{\sigma} - \hat{\sigma}_{est})^2]$.
# Estimateurs par maximum de vraissemblance
## Estimateurs par maximum de vraissemblance
Etant donné un modèle paramétrique et un jeu de données supposé avoir été généré par ce modèle, l'estimateur par *maximum de vraissemblance* ("maximum likelihood") associe à un paramètre la valeur qui rend le jeu de données observé le plus probable.
@ -33,26 +24,29 @@ Etant donné un modèle paramétrique et un jeu de données supposé avoir été
& \partial P(\mathbf{x}; \mu, \sigma^2) / \partial \sigma^2 = 0
\end{align*}
## Exemple d'une loi normale
### Exemple d'une loi normale
Supposons que des échantillons scalaires $X = x_1, x_2, \dots, x_n$ aient été générés selon une loi normale.
\begin{align*}
$$
\begin{aligned}
x_i &\sim \mathcal{N}(\mu, \sigma^2) \\
P(x_i ; \mu, \sigma) &= \frac{1}{\sqrt{2\pi\sigma^2}}exp\left( -\frac{1}{2\sigma^2}(x_i-\mu)^2\right)
\end{align*}
\end{aligned}
$$
\begin{align*}
$$
\begin{aligned}
& P(X ; \mu, \sigma) \\
= \{& \text{Hypothèse : les échantillons sont indépendants} \} \\
& \prod_i P(x_i ; \mu, \sigma) \\
= \phantom{\{}& \\
& (2\pi\sigma^2)^{-n/2} exp\left[ -\frac{1}{2\sigma^2} \sum_i (x_i-\mu)^2 \right]
\end{align*}
\end{aligned}
$$
Supposons que la moyenne $\hat{\mu}$ du modèle soit connue. Calculons alors l'estimateur par maximum de vraissemblance de la variance.
\begin{align*}
$$
\begin{aligned}
& \hat{\sigma}^2_{ML} \\
= \{& \text{Par définition d'un estimateur par maximum de vraissemblance.} \} \\
& \argmax_{\sigma^2} P(X ; \hat{\mu}, \sigma^2) \\
@ -70,11 +64,12 @@ Supposons que la moyenne $\hat{\mu}$ du modèle soit connue. Calculons alors l'e
& (n/2) s^{-2} \left[ \left( 1/n \sum_i (x_i - \hat{\mu})^2 \right) - s \right] = 0 \\
\Rightarrow \{& \text{Pour une variance finie, le second facteur doit s'annuler.} \} \\
& \hat{s}_{ML} = \hat{\sigma}^2_{ML} = 1/n \sum_i (x_i - \hat{\mu})^2
\end{align*}
\end{aligned}
$$
Quel est le biais de cet estimateur ?
\begin{align*}
$$
\begin{aligned}
& E\left[ \hat{\sigma}^2_{ML} \right] \\
= \{& \text{Voir dérivation ci-dessus.} \} \\
& E\left[ n^{-1} \sum_i (x_i - \hat{\mu})^2 \right] \\
@ -87,13 +82,14 @@ Quel est le biais de cet estimateur ?
& E[x_i^2] - 2 \hat{\mu}^2 + \hat{\mu}^2 \\
= \{& \hat{\sigma}^2 = E[x_i^2] - E[x_i]^2 = E[x_i^2] - \hat{\mu}^2 \} \\
& \hat{\sigma}^2
\end{align*}
\end{aligned}
$$
Cet estimateur est sans biais. On peut montrer qu'un estimateur non biaisé obtenu par l'approche du maximum de vraissemblance est également de variance minimale.
Considérons maintenant que la moyenne ne soit pas connue et calculons les estimateurs par maximum de vraissemblance de la moyenne et de la variance.
\begin{align*}
$$
\begin{aligned}
& \hat{\mu}_{ML} \\
= \{& \text{Par définition d'un estimateur par maximum de vraissemblance.} \} \\
& \argmax_{\mu} P(X ; \mu, s) \\
@ -108,15 +104,15 @@ Considérons maintenant que la moyenne ne soit pas connue et calculons les estim
& s^{-1} \left( \sum_i x_i \; - n \mu \right) = 0 \\
\Rightarrow \{& \text{Arithmétique} \} \\
& \hat{\mu}_{ML} = \frac{1}{n} \sum_i x_i
\end{align*}
\end{aligned}
$$
Quel est le biais de cet estimateur ?
$$E[\hat{\mu}_{ML}] = E[1/n \sum_i x_i] = 1/n \sum_i E[x_i] = 1/n \times n \times E[x_i] = \hat{\mu}$$
Il s'agit don d'un estimateur non biaisé.
\begin{align*}
$$
\begin{aligned}
& \hat{s}_{ML} \\
= \{& \text{Voir la précédente dérivation quand la moyenne était supposée connue.} \} \\
& n^{-1} \sum_i (x_i - \mu)^2 \\
@ -128,11 +124,12 @@ Il s'agit don d'un estimateur non biaisé.
& n^{-1} \sum_i x_i^2 \; - 2 \left( n^{-1} \sum_i x_i \right)^2 \; + \left( n^{-1} \sum_i x_i \right)^2 \\
= \{& \text{Arithmétique} \} \\
& n^{-1} \sum_i x_i^2 \; - \left( n^{-1} \sum_i x_i \right)^2
\end{align*}
\end{aligned}
$$
Quel est le biais de cet estimateur de la variance ?
\begin{align*}
$$
\begin{aligned}
& E[\hat{s}_{ML}] \\
= \{& \text{Voir dérivation ci-dessus.} \} \\
& E \left[ n^{-1} \sum_i x_i^2 \; - \left( n^{-1} \sum_i x_i \right)^2 \right] \\
@ -150,16 +147,17 @@ Quel est le biais de cet estimateur de la variance ?
& \hat{s} - n^{-1} \hat{s} \\
= \{& \text{Arithmétique} \} \\
& \frac{n-1}{n} \hat{s}
\end{align*}
\end{aligned}
$$
L'estimateur de la variance $\hat{s}_{ML}$ est maintenant biaisé. Nous pouvons construire un estimateur sans biais :
L'estimateur de la variance $\hat{s}_{ML}$ est maintenant biaisé. Nous pouvons construire un estimateur sans biais :
$$
s' \; = \; \left( \frac{n-1}{n} \right)^{-1} \hat{s}_{ML} \; = \; \frac{1}{n-1} \sum_i (x_i - \mu)^2
$$
$s'$ est alors sans biais mais n'est plus de variance minimale (bien que, parmi les estimateurs non biaisés, il soit de variance minimale).
# Analyse biais-variance pour la régression
## Analyse biais-variance pour la régression
Nous supposons que les données observées aient été générées par une fonction de la forme $y = f(\mathbf{x}) + \epsilon$ avec $\epsilon$ un bruit gaussien de moyenne nulle et de variance $\sigma^2$.
@ -168,8 +166,8 @@ A partir d'un jeu de données $\left\{ (\mathbf{x_i}, y_i) \right\}$, nous appre
Pour un nouveau point $\mathbf{x}^*$ qu'elle est l'espérance de l'erreur commise sur la prédiction de $y^* = f(\mathbf{x}^*) + \epsilon$, soit $E\left[\left(y^* - h(\mathbf{x}^*)\right)^2\right]$. Il s'agit de la moyenne de l'erreur sur l'ensemble infini de tous les jeux de données d'entraînement possibles.
Notons $\overline{x} = E[x]$, la valeur moyenne de x. Rappelons un résultat utile :
\begin{align*}
$$
\begin{aligned}
& E[(x - \overline{x})^2] \\
= \{& \text{Arithmétique} \} \\
& E[x^2 - 2 x \overline{x} + \overline{x}^2] \\
@ -179,11 +177,12 @@ Notons $\overline{x} = E[x]$, la valeur moyenne de x. Rappelons un résultat uti
& E[x^2] - 2 \overline{x}^2 + \overline{x}^2 \\
= \{& \text{Arithmétique} \} \\
& E[x^2] - \overline{x}^2
\end{align*}
\end{aligned}
$$
Nous décomposons l'espérance de l'erreur de prédiction ("Expected Prediction Error" ou EPE) d'un modèle de régression en biais, variance et bruit.
\begin{align*}
$$
\begin{aligned}
& E\left[ \left( h(\mathbf{x}^*) - y^* \right)^2 \right] \\
= \{& \text{Linéarité de E} \} \\
& E[h(\mathbf{x}^*)^2] - 2 E[h(\mathbf{x}^*)] E[y^*] + E[y^{*2}] \\
@ -194,7 +193,8 @@ Nous décomposons l'espérance de l'erreur de prédiction ("Expected Prediction
& E\left[ \left( h(\mathbf{x}^*) - \overline{h(\mathbf{x}^*)} \right)^2 \right] + \left( \overline{h(\mathbf{x}^*)} - f(\mathbf{x}^*) \right)^2 + \sigma^2 \\
= \{& \text{Introduction des définition de la variance, du biais et du bruit.} \} \\
& \text{Variance} \; + \; \text{Biais}^2 \; + \; \text{Bruit}^2
\end{align*}
\end{aligned}
$$
* La variance mesure la variation de la prédiction $h(\mathbf{x}^*)$ d'un jeu de données d'entraînement à l'autre.
* Le biais mesure l'erreur moyenne de $h(\mathbf{x}^*)$.

21
17_biais_variance_ridge.Rmd

@ -1,15 +1,6 @@
---
title: "17 Biais et variance des paramètres d'une régression ridge"
output:
bookdown::pdf_document2:
number_section: yes
includes:
in_header: preamble.tex
toc: false
classoption: fleqn
---
# Biais et variance des paramètres d'une régression ridge
# Biais des coefficients d'un modèle de régression ridge
## Biais des coefficients d'un modèle de régression ridge
Nous rappelons l'expression de l'estimation des coefficients d'un modèle de régression ridge en fonction de la décomposition en valeurs singulières de la matrice des données.
$$
@ -49,11 +40,11 @@ $$
Puisque $\sum \mathbf{v_j}\mathbf{v_j}^T=\mathbf{V}^T\mathbf{V}=\mathbf{I}_p$, dans le cas de la régression non régularisée (i.e., $\lambda=0$), $\hat{\boldsymbol\beta}$ est un estimateur sans biais (i.e., $E[\hat{\boldsymbol\beta}]=\boldsymbol\beta$). Mais pour la régression ridge, l'estimateur $\hat{\boldsymbol\beta}_\lambda$ est porteur de biais, ses composantes sont biaisées vers $0$. Elles le sont d'autant plus le long des directions de faible variance du jeu de données d'entraînement (c'est-à-dire les direction avec un faible $d_j$).
# Variance des coefficients d'un modèle de régression ridge
## Variance des coefficients d'un modèle de régression ridge
Avant de calculer la variance des coefficients d'un modèle de régression ridge, nous dérivons une forme utile de la variance des coefficients d'un modèle de régression linéaire non régularisé, toujours sous hypothèse de données générées selon un processus linéaire avec bruit gaussien de moyenne nulle (voir équation \ref{eq:processus-lineaire-bruit-gaussien}).
## Variance des coefficients d'un modèle de régression linéaire
### Variance des coefficients d'un modèle de régression linéaire
Nous utilisons la décomposition en valeurs singulières pour exprimer le modèle de l'équation \ref{eq:processus-lineaire-bruit-gaussien} en fonction d'une orthogonalisation de la matrice des données $\mathbf{X}$.
$$
@ -128,7 +119,7 @@ $$
\end{aligned}
$$
## Relation linéaire entre $\hat{\boldsymbol\beta}_\lambda$ et $\hat{\boldsymbol\beta}$
### Relation linéaire entre $\hat{\boldsymbol\beta}_\lambda$ et $\hat{\boldsymbol\beta}$
Nous allons découvrir une relation linéaire entre $\hat{\boldsymbol\beta}_\lambda$ et $\hat{\boldsymbol\beta}$. Elle nous permettra ensuite de calculer $Var\left(\hat{\boldsymbol\beta}_\lambda\right)$ sous hypothèse du modèle génératif (\ref{eq:processus-lineaire-bruit-gaussien}).
$$
@ -164,7 +155,7 @@ $$
\hat{\boldsymbol\beta}_\lambda = \mathbf{W} \hat{\boldsymbol\beta}
$$
## Variance des coefficients d'un modèle de régression linéaire régularisée
### Variance des coefficients d'un modèle de régression linéaire régularisée
Toujours sous hypothèse du modèle génératif (\ref{eq:processus-lineaire-bruit-gaussien}), nous pouvons maintenant calculer une expression de la variance des coefficients d'un modèle de régression linéaire régularisée en fonction de la décomposition en valeurs singulières de la matrice des données.
$$

63
18_kernel_ridge_regression.Rmd

@ -1,15 +1,12 @@
---
title: "18 Régression ridge à noyau"
output:
bookdown::pdf_document2:
number_section: yes
includes:
in_header: preamble.tex
toc: false
classoption: fleqn
---
# Régression ridge à noyau
# Noyau gaussien
```{r, include=FALSE}
source("01_intro.R", local = knitr::knit_global())
source("04_validation_croisee.R", local = knitr::knit_global())
source("18_kernel_ridge_regression.R", local = knitr::knit_global())
```
## Noyau gaussien
Soit un jeu de données :
$$
@ -51,15 +48,15 @@ xs <- seq(from=0, to=5, by=0.1)
fi <- function(ci,xi) function(x) ci * exp(-(x-xi)^2)
call_f <- function(f,...) f(...)
m <- t(matrix(unlist(lapply(Map(fi,ci,xi), call_f, xs)), nrow=npts, byrow=TRUE))
f <- rowSums(m)
plot(xs, f, type="l", lty="solid", ylim=range(cbind(f,m)))
g <- rowSums(m)
plot(xs, g, type="l", lty="solid", ylim=range(cbind(g,m)))
matplot(xs, m, type="l", lty="dotted", col=1, add=TRUE)
points(x=xi, y=rep(0,npts))
```
Cette approche est qualifiée de \emph{non paramétrique} car elle s'adapte aux données au lieu de chercher à adpater les paramètres d'une forme fonctionnelle fixée a priori comme nous le faisions auparavant dans le cas de la régression régularisée.
# Régression ridge à noyau
## Régression ridge à noyau
Nous montrons que l'approche non paramétrique introduite ci-dessus peut se présenter sous la forme d'une régression ridge après transformation des variables $\mathbf{x_i}$ par une fonction $\phi$.
$$
@ -78,7 +75,7 @@ $$
& \text{Identité de Woodbury} \\
= \{& \text{\url{https://en.wikipedia.org/wiki/Woodbury_matrix_identity}} \} \\
& (\mathbf{I} + \mathbf{U}\mathbf{V})^{-1}\mathbf{U} = \mathbf{U}(\mathbf{I} + \mathbf{V}\mathbf{U})^{-1} \\
= \{& \} \\
= \{& \text{Algèbre} \} \\
& \frac{\lambda}{\lambda}(\mathbf{I} + \mathbf{U}\mathbf{V})^{-1}\mathbf{U} = \frac{\lambda}{\lambda}\mathbf{U}(\mathbf{I} + \mathbf{V}\mathbf{U})^{-1} \\
= \{& \left(\lambda\mathbf{X}\right)^{-1} = \lambda^{-1}\mathbf{X}^{-1} \} \\
& (\lambda\mathbf{I} + \lambda\mathbf{U}\mathbf{V})^{-1}\lambda\mathbf{U} = \lambda\mathbf{U}(\lambda\mathbf{I} + \mathbf{V}\lambda\mathbf{U})^{-1} \\
@ -121,9 +118,9 @@ $$
$$
Ainsi, l'approche non paramétrique par juxtaposition des similarités d'une observation avec chaque observation du jeu d'entraînement correspond à une régression ridge après application de la transformation $\phi$ aux observations. Nous remarquons qu'il n'est pas nécessaire d'opérer explicitement cette transformation, seules sont nécessaires les valeurs des produits scalaires $\phi(\mathbf{x_i})^T\phi(\mathbf{x_j})$. En fait, l'application explicite de la transformation $\phi$ serait souvent impossible. Par exemple, dans le cas d'un noyau gaussien, $\phi$ projette vers un espace de dimension infinie...
# Calcul des paramètres de la régression ridge à noyau
## Calcul des paramètres de la régression ridge à noyau
Notons $\mathbf{G} = \mathbf{K} + \lambda \mathbf{I_n}$. Nous venons de montrer :
Notons $\mathbf{G} \triangleq \mathbf{K} + \lambda \mathbf{I_n}$. Nous venons de montrer :
\begin{equation}
\boldsymbol\alpha_\lambda = \mathbf{G}^{-1} \mathbf{y}
\label{eq:3}
@ -160,7 +157,7 @@ $$
$$
Ainsi, après avoir effectuée la décomposition en valeurs propres de $\mathbf{K}$ en $\mathcal{O}(n^3)$, le calcul de $\boldsymbol\alpha_\lambda$ pour une valeur $\lambda$ donnée se fait en $\mathcal{O}(n^2)$.
# LOOCV pour le choix de l'hyper-paramètre $\lambda$
## LOOCV pour le choix de l'hyper-paramètre $\lambda$
Pour la régression ridge, nous avons découvert une forme de la validation croisée un-contre-tous (LOOCV) qui peut être calculée à partir d'un seul calcul des paramètres sur tout le jeu d'entraînement :
$$
@ -226,6 +223,32 @@ $$
Nous pouvons donc calculer $\left(\mathbf{G}^{-1}\right)_{ii}$ en $\mathcal{O}(n)$ et $LOO_\lambda$ en $\mathcal{O}(n^2)$.
# Choix de l'hyper-paramètre $\sigma^2$
## Choix de l'hyper-paramètre $\sigma^2$
Il est possible de montrer qu'une fois les données standardisées (i.e., moyenne nulle et écart-type unité), la moyenne de la distance euclidienne entre deux points est égale à deux fois la dimension de l'espace des observations : $E\left[\|\mathbf{x_j}-\mathbf{x_i}\|^2\right] = 2d$. Ainsi, $\sigma^2=d$ est un choix raisonnable pour permettre au noyau gaussien de bien capturer les similarités entre points (i.e., les points les plus proches auront une similarité proche de $1$, tandis que les points les plus éloignés auront une similarité proche de $0$).
Il est possible de montrer qu'une fois les données standardisées (i.e., moyenne nulle et écart-type unité), la moyenne de la distance euclidienne entre deux points est égale à deux fois la dimension de l'espace des observations : $E\left[\|\mathbf{x_j}-\mathbf{x_i}\|^2\right] = 2d$. Ainsi, $\sigma^2=d$ est un choix raisonnable pour permettre au noyau gaussien de bien capturer les similarités entre points (i.e., les points les plus proches auront une similarité proche de $1$, tandis que les points les plus éloignés auront une similarité proche de $0$).
## Exemple sur un jeu de données synthétique
Nous reprenons le jeu de données synthétique utilisé depuis le premier module pour tester l'implémentation de la régression ridge à noyau gaussien.
```{r}
set.seed(1123)
n <- 100
data = gendat(n,0.2)
splitres <- splitdata(data,0.8)
entr <- splitres$entr
test <- splitres$test
km <- krr(entr$X,entr$Y)
yh <- predict(km,test$X)
plt(test,f)
points(test$X, yh, pch=4)
testmae <- mean(abs(yh-test$Y))
```
Ce modèle atteint une erreur absolue moyenne de `r testmae` sur le jeu de test.
## Annexe code source
```{r, code=readLines("18_kernel_ridge_regression.R"), eval=FALSE}
```

1
preamble.tex

@ -1,3 +1,4 @@
\usepackage[french]{babel}
\usepackage{amsmath}
\usepackage{url}
\DeclareMathOperator*{\argmax}{argmax}