Commençons par tester ce qu’on sait déjà faire, la régression logistique :

logreg <- glm(class ~ x1 + x2, family = binomial, data = spiral)
correct <- sum((fitted(logreg) > .5) + 1 == as.integer(spiral$class))
print(correct/250*100)

## [1] 56

140 points sur 250 sont bien classifiés, soit une précision de 56 % seulement.

Visuellement, ça donne ça :

C’est pas top ! cette performance n’est pas acceptable.

a. D’abord, l’architecture globale

Nous allons coder un réseau dans lequel les données progressent toujours en avant (feed-forward). l’information se déplace dans une seule direction, de la couche d’entrée à la couche de sortie en passant par les différentes couches cachées, dans l’ordre, comme dans le schéma précédent.

En termes de code, ça ressemble à ceci. On en profite pour faire en sorte que chaque neurone applique une sigmoïde (fonction logistique) juste avant de renvoyer sa sortie, on parle de fonction d’activation.

Nous allons écrire une fonction feedfwd() qui dans l’ordre : * Prend en entrée les variables explicatives (x), et deux vecteurs de poids (w1 et w2) * Multiplie x par w1 et lui applique la fonction d’activation * Le résultat de l’étape précédente subit la même opération avec le vecteur w2, puis est retourné. La sortie est un vecteur de probabilités.

# La fonction logistique 
sig <- function(x) 1 / (1 + exp(-x))

# L'opérateur %*% sert à la multiplication de matrices
feedfwd <- function(x, w1, w2) {
  z1 <- cbind(1, x) %*% w1
  h <- sig(z1)
  z2 <- cbind(1, h) %*% w2
  list(output = sig(z2), h = h)
}

b. Ensuite, rajoutons un peu d’intelligence avec la rétro-propagation du gradient

La rétro-propagation du gradient est le cœur de la machinerie permettant d’entraîner un réseau de neurones. C’est elle qui permet progressivement de trouver les poids \(W\) optimaux qui font que notre modèle est “bon”. La rétro-propagation permet au modèle d’apprendre de ses erreurs de manière itérative et de s’améliorer au fur et à mesure qu’il ingère les données. Les poids qui contribuent à engendrer une erreur importante se verront modifiés de manière plus significative que les poids qui ont engendré une erreur marginale, dans la direction de la minimisation de l’erreur globale du réseau. L’algorithme du gradient a pour but de converger de manière itérative vers une configuration optimale des poids, permettant le meilleur ajustement de modèle possible.

Quelques points à ce sujet :

• Tout algorithme d’apprentissage supervisé a pour objectif d’apprendre une structure de données à partir d’une base d’apprentissage. Pour ce faire, il a besoin de trouver la combinaison de poids qui permet la meilleure approximation des données d’apprentissage.
• Pour atteindre ce meilleur ajustement, on optimise (= minimise la valeur) d’une fonction de coût (Objectif/Loss/Cost). Cette fonction de coût mesure l’écart entre les vrais données \(y\) et les prédictions du modèle \(\hat{y}\), c’est à dire l’erreur. Nous sommes bien face à un problème d’optimisation avec une fonction de coût à minimiser.
• Ce problème d’optimisation, en apparence complexe, se résout en fait assez efficacement grâce à l’algorithme de la descente du gradient. Le gradient est une généralisation de la notion de dérivée, dans le cas à plusieurs variables (vous trouverez ici une explication imagée et intuitive de l’algo). Dans l’espace à deux dimensions (2 variables), voici une illustration de l’algorithme (source) :
  • On a ici deux coefficients de régression m et b à chercher de manière à ce que notre droite apprenne les données (encart de droite).
  • Pour ce faire, on cherche le minimum de Error (encart de gauche) :

Hop :

backprop <- function(x, y, y_hat, w1, w2, h, learn_rate) {
  dw2 <- t(cbind(1, h)) %*% (y_hat - y)
  dh  <- (y_hat - y) %*% t(w2[-1, , drop = FALSE])
  dw1 <- t(cbind(1, x)) %*% (h * (1 - h) * dh)
  
  w1 <- w1 - learn_rate * dw1
  w2 <- w2 - learn_rate * dw2
  
  list(w1 = w1, w2 = w2)
}

c. Rassemblons tout ça et testons !

On définit donc une fonction train() permettant d’entraîner notre réseau de neurones, avec une architecture feed-forward et utilisant la rétro-propagation pour trouver les meilleurs paramètres. Les paramètres sont :

• x : matrice comprenant les données d’apprentissage (les variables x1 et x2)
• y : vecteur contenant la variable réponse (label : “Rouge” ou “Bleu”)
• hidden : nombre de nœuds cachés de la couche cachée
• learn_rate : learning rate ou taux d’apprentissage. Ce paramètre contrôle le pas avec lequel s’effectue la descente du gradient. Plus il est grand, et plus la descente se fait à grands pas et donc rapidement (en théorie), mais on prend le risque de “sauter” le minimum de la fonction. Prenons donc une petite valeur (par défaut 0.01)
• iterations : nombre d’itérations (une itération correspond à un passage des données dans le réseau suivi d’une mise à jour des coefficients par rétro-propagation)

Les poids w sont initialisés avec des valeurs aléatoires, donc pas besoin de les fournir à la fonction train() et on ré-utilise la fonction feedfwd() créée précédemment

train <- function(x, y, hidden = 5, learn_rate = 1e-2, iterations = 1e4) {
  d <- ncol(x) + 1
  w1 <- matrix(rnorm(d * hidden), d, hidden) 
  w2 <- as.matrix(rnorm(hidden + 1))
  for (i in 1:iterations) {
    ff <- feedfwd(x, w1, w2)
    bp <- backprop(x, y,
                        y_hat = ff$output,
                        w1, w2,
                        h = ff$h,
                        learn_rate = learn_rate)
    w1 <- bp$w1; w2 <- bp$w2
  }
  list(output = ff$output, w1 = w1, w2 = w2)
}

Essayons avec le jeu de données en spirales avec une couche cachée composée de 5 nœuds et en faisant 100 000 itérations :

x <- data.matrix(spiral[, c('x1', 'x2')])
y <- spiral$class == 'Bleu'
nnet5 <- train(x, y, hidden = 5, iterations = 1e5)

En termes de performance de prédiction, on est à 98 %, ce qui est vraiment pas mal avec des données complexes à saisir et un entrainement très rapide (10s sur un processeur M1)

Graphiquement, il est intéressant d’observer la forme de la frontière de décision, non linéaire mais assez anguleuse quand même. On commence à percevoir la puissance d’apprentissage des réseaux de neurones :

On entraîne à présent un réseau avec une couche cachée plus grande, constituée de 30 nœuds :

nnet30 <- train(x, y, hidden = 30, iterations = 1e5)

ff_grid <- feedfwd(x = data.matrix(grid[, c('x1', 'x2')]),
                       w1 = nnet30$w1,
                       w2 = nnet30$w2)
grid$class <- factor((ff_grid$output > .5) * 1,
                     labels = levels(spiral$class))

La frontière de décision est totalement arrondie, parfaitement adaptée à nos données. La performance est de … 100% !!!

Testons à présent le cas avec un seul neurone dans la couche cachée :

nnet1 <- train(x, y, hidden = 1, iterations = 1e5)

Nous voilà de retour au cas de la régression linéaire, une frontière de décision linéaire et une performance médiocre :

ff_grid <- feedfwd(x = data.matrix(grid[, c('x1', 'x2')]),
                       w1 = nnet1$w1,
                       w2 = nnet1$w2)
grid$class <- factor((ff_grid$output > .5) * 1,
                     labels = levels(spiral$class))