Le document suivant a été réalisé avec le logiciel libre https://www.r-project.org/. Les informations sur les versions précises du logiciel et des packages ainsi que sur la configuration de la machine (en particulier, ce document a été réalisé sous Linux et la prise en charge des caractères accentués est réalisée par l’encodage UTF-8 ; pour compiler correctement cette vignette sous Windows, il sera probablement nécessaire d’utiliser le menu « Fichier / Ré-ouvrir avec l’encodage … ») au moment de la compilation sont fournies à la fin de ce document (section « Information sur la session R »).

Le but de ce document est d’illustrer, sur des cas très simples, l’utilisation et le comportement des perceptrons à une couche cachée. Les études de cas sont réalisées avec les packages nnet (apprentissage du réseaux de neurones) et e1071 (calibration des divers hyper-paramètres de la méthode). En outre, le package mlbench est utilisé pour fournir un exemple de données réelles (dans la deuxième partie de ce document « Un exemple sur des données réelles ») et le package devtools, non indispensable pour les applications, est simplement utilisé pour fournir des informations sur le logiciel, les packages et l’environnement de compilation utilisés.

Les packages peuvent être installés avec la commande

install.packages(c("nnet", "e1701", "mlbench", "devtools"))

et ils sont chargés dans l’environnement de travail de R avec

library(devtools)
library(nnet)
library(e1071)
library(mlbench)

Enfin, le présent document (format HTML) ainsi que sa source (format RMarkdown, compilable en utilisant RStudio et les packages evaluate, formatR, highr, markdown, yaml, htmltools, knitr et rmarkdown ; ces packages s’installent automatiquement à la première compilation d’un fichier au format RMarkdown sous RStudio) sont disponible sur mon site web.

Données simulées

La première partie de ce document illustre le comportement du perceptron dans une tâche simple qui est l’estimation d’une fonction de R dans R à partir d’observations bruitées.

Les données sont générées par les commandes suivantes (dans lesquelles la fonction set.seed est utilisée pour initialiser une graine aléatoire et assurer des résultats reproductibles) :

set.seed(1528)
x <- runif(15, 0, 1)
y <- sin(x*2*pi) + rnorm(15, 0, 0.3)

Les données générées sont représentées dans la figure ci-dessous avec la fonction sous-jacente à estimer, \(x \rightarrow \sin(2\pi x)\) :

t <- seq(0, 1, length=1000)
yt <- sin(t*2*pi)
plot(x, y, pch = "+", lwd = 2)
lines(t, yt, col = "darkred")

Influence de \(Q\)

Dans un premier temps, nous illustrons l’influence de l’hyper-paramètre \(Q\) (le nombre de neurones sur la couche cachée). L’apprentissage du perceptron à une couche cachée est réalisée avec la fonction nnet sur laquelle on peut obtenir de l’aide en utilisant la commande

help(nnet)

Dans ce qui suit, pour les observations générées précédemment et pour chaque valeur de \(Q\) variant entre 1 et 10, 10 apprentissages sont réalisés. Le perceptron qui minimise l’erreur quadratique empirique sur les données d’apprentissage est finalement conservé (un perceptron pour chaque valeur de \(Q\)) :

set.seed(1513)
Qlist <- 1:10
res <- list()
iter <- 10
for (k in Qlist) {
  best_error <- 10^5
  for (ind in 1:iter) {
    cur_nnet <- nnet(x, y, size = k, trace = FALSE, linout = TRUE)
    if (mean(cur_nnet$residuals^2 < best_error))
      best_nnet <- cur_nnet
  }
  res[[k]] <- best_nnet
}

La fonction estimée par le perceptron est ensuite représentée en comparaison de la fonction cible. Pour cela, on utilise la fonction predict qui permet d’obtenir la prédiction du perceptron pour 1000 valeurs de \(x\) uniformément réparties dans \([0,1]\). L’erreur moyenne réalisée par le perceptron sur ces 100 valeurs est conservées dans la variable gen_err : cette erreur moyenne sera comparée à l’erreur moyenne d’apprentissage (conservée dans la variable train_err) :

par(mfrow = c(2, 5))
train_err <- unlist(lapply(res, function(anet) anet$value/length(x)))
gen_err <- vector(length = 10)
for (k in Qlist) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[k]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = paste0("Q = ", k),
       ylim = range(c(y, yt, cur_pred)), cex=2, cex.main=2, axes=FALSE)
  box()
  lines(t, yt, col = "darkred", lwd=2, lty=3)
  lines(t, cur_pred, col = "black", lwd=2)
  gen_err[k] <- mean((cur_pred - yt)^2)
}

Le graphique ci-dessous présente justement une comparaison de l’évolution de l’erreur empirique d’apprentissage (train_err) avec l’erreur commise sur la grille régulière de 1000 points sur \([0,1]\) (gen_err) qui sera appelée « erreur de test » et correspond à une estimation de la quantité \[ \mathbb{E}\left[\left(\hat{f}_Q(X) - f(X) \right)^2\right] \] pour \(X \sim \mathcal{U}[0,1]\) et \(f(x) = \sin(2\pi x)\) :

plot(Qlist, train_err, pch = "+", col = "darkred", type = "b", 
     ylim = c(0, max(c(train_err, gen_err))), xlab = "Q", ylab = "erreur", 
     cex.axis = 1.5, cex.lab = 1.5)
lines(Qlist, gen_err, col = "black", lty=2)
points(Qlist, gen_err, pch = "+", col = "black")
legend(0.7, 0.25, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))

On y observe que l’erreur d’apprentissage a tendance à décroître lorsque \(Q\) augmente alors que l’erreur de test a une phase de décroissance suivi (grossièrement) d’une phase de croissance. Le minimum de l’erreur de test est atteint pour \(Q=\) 3 avec une erreur égale à :

min(gen_err)
## [1] 0.04503842

Calibration de \(Q\)

En l’absence d’un grand échantillon de test pour calibrer la valeur de \(Q\) par calcul d’une erreur de tests, diverses stratégies classiques (validation simple, validation croisée, bootstrap…) peuvent être mises en œuvre pour sélectionner une valeur pertinente pour \(Q\). Ces stratégies peuvent être mise en œuvre grâce à la fonction tune.nnet du package e1071. Les diverses options possibles pour la calibration (type de calibration - par validation simple, croisée, bootstrap … - ou paramètre de la méthode utilisée comme le nombre de groupes dans la validation croisée) sont controlées via la fonction tune.control dont le résultat est passé en argument à la fonction tune.nnet.

Dans la commande ci-dessous, on illustre la callibration par validation croisée (10 groupes) de la valeur de \(Q\) pour les données précédentes. Ici, l’option nrepeat sert à fixer le nombre d’initialisations différentes (et donc de perceptrons différents) qui sont effectuées pour produire, pour chaque groupe de la validation croisée, une estimation différente de l’erreur :

set.seed(1634)
res_tune <- tune.nnet(x, y, tune.control(nrepeat = iter), size = Qlist, 
                      linout = TRUE)

Le graphique ci-dessous illustre l’évolution de l’erreur de validation croisée obtenue pour les différentes valeurs de \(Q\) :

plot(Qlist, res_tune$performances$error, main = "", xlab = "Q", ylab = "erreur",
     type = "b")

L’erreur de validation croisée est donc minimale pour \(Q=\) 1 et on peut obtenir les informations additionnelles suivantes sur sa valeur et sur les caractéristiques du perceptron obtenu à partir de ce paramètre optimal :

res_tune$best.parameters
##   size
## 1    1
res_tune$best.performance
## [1] 0.1263338
res_tune$best.model
## a 1-1-1 network with 4 weights
## options were - linear output units

Régularisation

Cette partie illustre l’utilisation d’une pénalité dans l’apprentissage du perceptron. De manière plus précise, la fonction objectif minimisée lors de l’apprentissage est : \[ \hat{\mathcal{R}}_n(w) + \lambda \|w\|^2 \] où \(\lambda\) est le paramètre de régularisation. Dans les lignes de commande ci-dessous, on fait varier ce paramètre selon une échelle exponentielle pour un perceptron à une couche cachée dont le nombre de neurones est très supérieur à l’optimal (on a choisi \(Q=10\)) :

set.seed(1546)
k <- 10
reg_list <- c(10^(- (1:9)), 0)
res <- list()
iter <- 10
for (ind in seq_along(reg_list)) {
  best_error <- 10^5
  for (repet in 1:iter) {
    cur_nnet <- nnet(x, y, size = k, trace = FALSE, linout = TRUE, 
                     decay = reg_list[ind])
    if (mean(cur_nnet$residuals^2 < best_error))
      best_nnet <- cur_nnet
  }
  res[[ind]] <- best_nnet
}

Les résultats obtenus, comme dans la partie « Calibration de \(Q\) » sont illustrées par le graphique suivant dans lequel on constate qu’un fort paramètre de régularisation conduit à une estimation trop grossière de la fonction cible alors qu’un faible paramètre de régularisation produit des résultats similaires à ceux obtenus sans régularisation, à savoir un sur-apprentissage lors de l’apprentissage :

par(mfrow = c(2,5))
train_err <- unlist(lapply(res, function(anet) anet$value/length(x)))
gen_err <- vector(length = 10)
print_reg <- c(reg_list[-10], "sans régularisation")
for (ind in seq_along(reg_list)) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[ind]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = print_reg[ind], 
       ylim = range(c(y, yt, cur_pred)), cex.main = 1, axes = FALSE)
  box()
  lines(t, yt, col = "darkred", lty=3, lwd=2)
  lines(t, cur_pred, col = "black", lwd=2)
  gen_err[ind] <- mean((cur_pred - yt)^2)
}

Le graphique ci-dessous illustre précisément l’évolution des erreurs d’apprentissage et de test selon la valeur de \(-\log_{10}(\lambda)\) (pour \(\lambda > 0\)) et montre qu’un minimum est atteint pour \(\lambda=\) 0.001 :

plot(-log10(reg_list[-10]), train_err[-10], pch = "+", col = "darkred",
     type = "b", ylim = c(0, max(c(train_err, gen_err))), 
     xlab = expression(-log[10](lambda)), ylab = "erreur", cex.axis = 1.5,
     cex.lab = 1.5)
lines(-log10(reg_list[-10]), gen_err[-10], col = "black", lty=2)
points(-log10(reg_list[-10]), gen_err[-10], pch = "+", col = "black")
legend(0.7, 0.3, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))

La valeur de l’erreur de test obtenue pour cet optimum est donnée ci-dessous. On constatera, en particulier, qu’elle est inférieur à l’erreur de test obtenue pour la valeur optimale de \(Q\) sans régularisation :

min(gen_err)
## [1] 0.03393102

Calibration du paramètre de régularisation

De la même manière que la calibration du nombre de neurones sur la couche cachée a été réalisée, calibrez le paramètre de régularisation \(\lambda\) par validation croisée en utilisant tune.nnet.

Si vous initialisez la graine aléatoire avec set.seed(1612), vous devriez obtenir les résultats suivants :

res_tune$best.parameters
##   size decay
## 3   10 0.001
res_tune$best.performance
## [1] 0.09905401
res_tune$best.model
## a 1-10-1 network with 31 weights
## options were - linear output units  decay=0.001

Un exemple sur des données réelles

Enfin, cette dernière partie vous propose de mettre en œuvre ce que vous avez appris sur les données BostonHousing du package mlbench. Ces données peuvent être chargées avec les commandes :

data("BostonHousing")

et on trouve des informations sur leur provenance et les questions qu’elles soulèvent dans l’aide :

?BostonHousing

En particulier, on cherchera à prédire le revenu médian des habitants des foyers d’un quartier (medv) à partir de caractéristiques socio-économiques de ce quartier. Pour simplifier le problème, nous retirons la variable explicative non quantitative (rad) et effectuons une standardisation des données :

df <- BostonHousing
df$chas <- NULL
df <- as.data.frame(scale(df))

Utilisez la fonction tune.nnet pour calibrer simultanément le nombre de neurones sur la couche cachée parmi \(\{5,\, 10,\, 15,\, 20\}\) et le paramètre de régularisation parmi \(\{1,\, 10^{-1},\, 10^{-2},\, 10^{-3}\}\). En fixant la graine aléatoire à la valeur set.seed(1409), vous obtiendrez les résultats suivants (pour une valeur de maxit égale à 1000 et 5 répétitions de l’apprentissage : attention ! l’apprentissage est long…) :

plot(res_tune, main = "", color.palette = topo.colors, transform.y = log10)

res_tune$best.parameters
##   size decay
## 9   20   0.1
res_tune$best.model
## a 12-20-1 network with 281 weights
## inputs: crim zn indus nox rm age dis rad tax ptratio b lstat 
## output(s): medv 
## options were - linear output units  decay=0.1
res_tune$best.model$convergence
## [1] 1

La qualité de l’ajustement peut être visualisée dans le graphique suivant :

plot(df$medv, res_tune$best.model$fitted.values, pch = 19, col = "darkred",
     cex = 0.7, xlab = "valeurs réelles", ylab = "valeurs prédites")
abline(0, 1)

Information sur la session R

session_info()
## Session info --------------------------------------------------------------
##  setting  value                       
##  version  R version 3.3.1 (2016-06-21)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language en_US                       
##  collate  en_US.UTF-8                 
##  tz       Europe/Paris                
##  date     2016-09-14
## Packages ------------------------------------------------------------------
##  package   * version date       source        
##  class       7.3-14  2015-08-30 CRAN (R 3.2.5)
##  codetools   0.2-14  2015-07-15 CRAN (R 3.3.1)
##  devtools  * 1.12.0  2016-06-24 CRAN (R 3.3.1)
##  digest      0.6.10  2016-08-02 CRAN (R 3.3.1)
##  e1071     * 1.6-7   2015-08-05 CRAN (R 3.3.1)
##  evaluate    0.9     2016-04-29 CRAN (R 3.3.1)
##  formatR     1.4     2016-05-09 CRAN (R 3.3.1)
##  htmltools   0.3.5   2016-03-21 CRAN (R 3.3.1)
##  knitr       1.14    2016-08-13 CRAN (R 3.3.1)
##  magrittr    1.5     2014-11-22 CRAN (R 3.3.1)
##  memoise     1.0.0   2016-01-29 CRAN (R 3.3.1)
##  mlbench   * 2.1-1   2012-07-10 CRAN (R 3.3.1)
##  nnet      * 7.3-12  2016-02-02 CRAN (R 3.3.1)
##  Rcpp        0.12.7  2016-09-05 CRAN (R 3.3.1)
##  rmarkdown   1.0     2016-07-08 CRAN (R 3.3.1)
##  stringi     1.1.1   2016-05-27 CRAN (R 3.3.1)
##  stringr     1.1.0   2016-08-19 CRAN (R 3.3.1)
##  withr       1.0.2   2016-06-20 CRAN (R 3.3.1)
##  yaml        2.1.13  2014-06-12 CRAN (R 3.3.1)