Réseaux de neurones et apprentissage profond

Aux chapitres 2 et 3, nous avons construit des modèles linéaires pour la régression et la classification. Au chapitre 4, nous avons vu comment enrichir ces modèles en transformant les entrées par une fonction $\boldsymbol{\phi}$ fixée à l’avance. Ce chapitre franchit une étape supplémentaire: au lieu de choisir $\boldsymbol{\phi}$ manuellement, nous allons l’apprendre à partir des données. Cette idée conduit aux réseaux de neurones.

Dans ce chapitre, nous rappelons d’abord le cadre probabiliste qui unifie régression et classification, puis nous montrons comment le perceptron simple atteint une limite structurelle (illustrée par le problème XOR). Nous présentons ensuite l’anatomie d’un réseau multicouche (couches, activations, architecture). La section sur la dérivation automatique est plus technique: elle développe la règle de la chaîne, les produits jacobien-vecteur, et l’algorithme de rétropropagation, puis montre comment les bibliothèques modernes implémentent ces idées. Vous pouvez survoler les détails en première lecture et retenir le mécanisme général. La section d’implémentation propose un MLP complet en NumPy, et le chapitre se termine par une mise en perspective montrant comment les MLP sont utilisés en pratique pour la régression et la classification. Les algorithmes d’optimisation, la stabilisation de l’entraînement et la régularisation sont couverts au chapitre suivant.

Le cadre unifié: prédire les paramètres d’une distribution¶

Régression et classification comme maximum de vraisemblance¶

Revenons au cadre probabiliste des chapitres 2 et 5. Dans tous les modèles que nous avons vus, le problème d’apprentissage supervisé prend la même forme: étant donné une entrée $\mathbf{x}$, nous voulons prédire les paramètres d’une distribution conditionnelle $p(y | \mathbf{x}; \boldsymbol{\theta})$, puis trouver $\boldsymbol{\theta}$ par maximum de vraisemblance.

En régression, nous avons supposé un bruit gaussien:

Le modèle prédit la moyenne $\mu(\mathbf{x})$ de la distribution. La log-vraisemblance négative donne, à une constante près, la perte des moindres carrés:

En classification binaire, nous avons supposé une distribution de Bernoulli:

Le modèle prédit la probabilité $\mu(\mathbf{x}) = p(y = 1 | \mathbf{x})$. La log-vraisemblance négative donne l’entropie croisée binaire. Pour la classification multiclasse, la distribution catégorielle et la fonction softmax jouent le même rôle, avec l’entropie croisée catégorielle comme perte.

Dans chaque cas, une fonction $\mu(\mathbf{x})$ prend une entrée et produit les paramètres de la distribution de sortie. Toute la question est: quelle forme donner à cette fonction?

Des modèles linéaires aux caractéristiques apprises¶

Jusqu’ici, nos modèles ont été linéaires dans les entrées. Pour la régression:

Pour la classification binaire, la probabilité passe par une sigmoïde, mais la pré-activation reste linéaire:

Au chapitre 4, nous avons étendu cette approche avec l’expansion de caractéristiques. Au lieu d’utiliser $\mathbf{x}$ directement, nous le transformons par une fonction $\boldsymbol{\phi}: \mathbb{R}^d \to \mathbb{R}^D$ choisie à l’avance (polynômes, fonctions trigonométriques, bases radiales, etc.):

Le modèle reste linéaire dans les paramètres $\boldsymbol{\theta}$, ce qui facilite l’optimisation, mais il capture des relations non linéaires en $\mathbf{x}$ grâce au choix de $\boldsymbol{\phi}$.

Cette approche a une limite importante: le choix de $\boldsymbol{\phi}$ repose entièrement sur l’expertise du praticien. Pour des données tabulaires simples, cela peut fonctionner. Mais pour des images, du texte ou de l’audio, concevoir manuellement les bonnes caractéristiques est très difficile, et souvent le facteur limitant de la performance.

Les réseaux de neurones paramètrent $\boldsymbol{\phi}$ et l’apprennent à partir des données. Au lieu d’écrire $\boldsymbol{\theta}^\top \boldsymbol{\phi}(\mathbf{x})$ avec $\boldsymbol{\phi}$ fixé, nous écrivons:

où $\boldsymbol{\phi}(\cdot; \boldsymbol{\theta}_\phi)$ est elle-même une fonction paramétrique. Les paramètres $\boldsymbol{\theta}_\phi$ contrôlent la transformation des entrées (les “caractéristiques apprises”), tandis que $\mathbf{w}$ sont les poids de la couche de sortie. On optimise les deux simultanément: le modèle apprend la représentation et le prédicteur en même temps.

Cela soulève deux questions: quelle forme donner à $\boldsymbol{\phi}(\cdot; \boldsymbol{\theta}_\phi)$, et comment optimiser l’ensemble des paramètres? Le reste de ce chapitre répond à ces deux questions.

Le perceptron: aux origines des réseaux de neurones¶

La régression logistique et les réseaux de neurones partagent un ancêtre commun: le perceptron, proposé par Rosenblatt en 1958 Rosenblatt (1958). Comprendre ce modèle et sa limitation éclaire pourquoi les réseaux multicouches ont été inventés, et pourquoi leur développement a pris plusieurs décennies.

Un modèle inspiré du neurone biologique¶

En 1943, McCulloch et Pitts McCulloch & Pitts (1943) ont formalisé le comportement du neurone biologique: une unité qui reçoit des signaux pondérés et s’active si leur somme dépasse un seuil. Le modèle se résume à:

Ce neurone calcule une combinaison linéaire des entrées, puis prend une décision binaire: actif ($\hat{y} = 1$) ou inactif ($\hat{y} = 0$). Rosenblatt y a ajouté un algorithme pour ajuster les poids $\boldsymbol{\theta}$ à partir d’exemples étiquetés. L’enthousiasme de l’époque était considérable: des démonstrateurs matériels ont été construits, et la presse grand public annonçait une machine capable «d’apprendre à reconnaître».

Le point de départ est la neuroscience computationnelle, pas l’optimisation. McCulloch et Pitts voulaient formaliser le comportement des neurones corticaux; Rosenblatt s’inspirait de la vision artificielle et des réseaux nerveux. C’est cette origine qui distingue la trajectoire intellectuelle du perceptron de celle de la régression logistique, même si les deux modèles aboutissent à une structure mathématique très proche.

Lien avec la régression logistique¶

Les deux modèles calculent la même pré-activation linéaire $z = \boldsymbol{\theta}^\top \mathbf{x}$, mais diffèrent dans la façon dont ils l’interprètent:

La frontière de décision est dans les deux cas le même hyperplan $\{\mathbf{x} : \boldsymbol{\theta}^\top \mathbf{x} = 0\}$. La sigmoïde peut être vue comme une version lisse et probabiliste de la fonction échelon: au lieu de trancher brusquement, elle exprime l’incertitude via une probabilité.

Source

import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# --- Panneau gauche: activations ---
z = np.linspace(-5, 5, 500)
sigmoid = 1 / (1 + np.exp(-z))

ax = axes[0]
ax.plot(z, sigmoid, '#1f77b4', lw=2.5, label=r'Sigmoïde $\sigma(z)$')
ax.plot(z[z < 0],  np.zeros(np.sum(z < 0)),  '#d62728', lw=2.5)
ax.plot(z[z >= 0], np.ones(np.sum(z >= 0)),  '#d62728', lw=2.5,
        label=r'Échelon $\mathbf{1}[z \geq 0]$')
ax.scatter([0], [1], color='#d62728', s=55, zorder=5)
ax.scatter([0], [0], color='white', edgecolors='#d62728', linewidths=2, s=55, zorder=5)
ax.axvline(0, color='gray', lw=1, linestyle=':', alpha=0.6)
ax.set_xlabel(r'Pré-activation $z = \boldsymbol{\theta}^\top \mathbf{x}$')
ax.set_ylabel('Sortie')
ax.set_title('Activations comparées')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5); ax.set_ylim(-0.1, 1.2)
ax.annotate(r'Frontière: $z=0$', xy=(0, 0.5), xytext=(1.8, 0.25),
            fontsize=9, color='#555555',
            arrowprops=dict(arrowstyle='->', color='gray', lw=1.2))

# --- Panneau droit: frontière de décision 2D ---
rng = np.random.default_rng(0)
n = 40
X_pos = rng.multivariate_normal([ 1.5,  0.8], [[0.4, 0], [0, 0.4]], n)
X_neg = rng.multivariate_normal([-1.5, -0.8], [[0.4, 0], [0, 0.4]], n)

ax = axes[1]
xx, yy = np.meshgrid(np.linspace(-3.5, 3.5, 200), np.linspace(-2.5, 2.5, 200))
Z = xx + yy   # theta = [1, 1], frontière: x1 + x2 = 0
ax.contourf(xx, yy, Z, levels=[-100, 0, 100],
            colors=['#fdd8d8', '#d8e8fd'], alpha=0.45)
ax.contour(xx, yy, Z, levels=[0], colors='k', linewidths=2)
ax.scatter(X_pos[:, 0], X_pos[:, 1], color='#1f77b4', marker='o', s=35,
           alpha=0.85, label='Classe 1', edgecolors='white', linewidths=0.5)
ax.scatter(X_neg[:, 0], X_neg[:, 1], color='#d62728', marker='s', s=35,
           alpha=0.85, label='Classe 0', edgecolors='white', linewidths=0.5)
ax.text( 1.8, -2.0, 'Rég. log.: prob. $\\to 1$\nPerceptron: classe 1',
         fontsize=7.5, color='#1f77b4', ha='center',
         bbox=dict(boxstyle='round,pad=0.3', fc='white', ec='#1f77b4', alpha=0.8))
ax.text(-1.8,  1.8, 'Rég. log.: prob. $\\to 0$\nPerceptron: classe 0',
         fontsize=7.5, color='#d62728', ha='center',
         bbox=dict(boxstyle='round,pad=0.3', fc='white', ec='#d62728', alpha=0.8))
ax.set_xlabel(r'$x_1$'); ax.set_ylabel(r'$x_2$')
ax.set_title(r'Même frontière: $\boldsymbol{\theta}^\top \mathbf{x} = 0$')
ax.legend(fontsize=9, loc='upper left')
ax.grid(True, alpha=0.3)

plt.suptitle('Régression logistique et perceptron: deux lectures du même hyperplan', fontsize=10)
plt.tight_layout()

La règle d’apprentissage et la difficulté de l’optimisation¶

La régression logistique minimise l’entropie croisée, une fonction différentiable, ce qui autorise la descente de gradient. Le perceptron minimise la perte perceptron (avec des étiquettes $y_i \in \{-1, +1\}$):

Cette perte est nulle pour les exemples bien classés et pénalise les erreurs proportionnellement à leur amplitude. Elle n’est toutefois pas différentiable au point exact où $\boldsymbol{\theta}^\top \mathbf{x}_i = 0$. On utilise alors le sous-gradient, qui généralise le gradient aux fonctions non différentiables:

Cela donne la règle de mise à jour classique: pour chaque exemple mal classé, corriger les poids dans la direction de cet exemple,

La convergence est plus délicate qu’avec la descente de gradient sur une fonction convexe et différentiable. Si les données sont linéairement séparables, le théorème de convergence du perceptron Novikoff (1962) garantit que l’algorithme s’arrête en un nombre fini de mises à jour, borné par $(R/\gamma)^2$ où $R$ est le rayon des données et $\gamma$ la marge de séparation. Si les données ne sont pas linéairement séparables, l’algorithme peut cycler sans jamais converger.

Une limite structurelle¶

Toutes ces variantes (perceptron, régression logistique, moindres carrés) partagent la même contrainte: leur frontière de décision est un hyperplan. Quelle que soit la façon dont on choisit ou entraîne les poids $\boldsymbol{\theta}$, on ne peut séparer que des classes linéairement séparables.

C’est précisément ce que Minsky et Papert ont formalisé en 1969 Minsky & Papert (1969), en montrant que certaines fonctions booléennes sont impossibles à apprendre pour un perceptron simple. L’exemple canonique est la fonction XOR (ou exclusif):

Les points de classe 0 sont disposés en diagonale, $(0,0)$ et $(1,1)$, et ceux de classe 1 sur l’autre, $(0,1)$ et $(1,0)$. Aucune droite ne peut séparer ces deux groupes. Mais Minsky et Papert montraient aussi la solution: empiler deux couches de perceptrons suffit, car la première couche peut transformer l’espace de sorte que les classes deviennent linéairement séparables. Nous verrons dans la section d’implémentation qu’un petit MLP résout XOR sans difficulté.

Leur analyse a contribué à un ralentissement de la recherche sur les réseaux de neurones pendant plusieurs années, jusqu’à ce que les avancées en optimisation et en calcul redonnent vie au domaine. La section suivante formalise l’architecture multicouche qui résout cette limitation.

Anatomie d’un réseau de neurones¶

Un neurone: transformation affine et non-linéarité¶

Un réseau de neurones est construit à partir d’une opération élémentaire: une transformation affine suivie d’une fonction non linéaire. Pour une entrée $\mathbf{x} \in \mathbb{R}^d$:

où $\mathbf{w} \in \mathbb{R}^d$ est un vecteur de poids, $b \in \mathbb{R}$ est un biais, et $\varphi: \mathbb{R} \to \mathbb{R}$ est une fonction d’activation non linéaire. La quantité $a = \mathbf{w}^\top \mathbf{x} + b$ est la pré-activation et $h$ est l’activation du neurone.

Une couche de $m$ neurones applique cette opération en parallèle, ce qui s’écrit sous forme matricielle:

où $W \in \mathbb{R}^{m \times d}$ est la matrice de poids, $\mathbf{b} \in \mathbb{R}^m$ le vecteur de biais, et $\varphi$ est appliquée élément par élément.

On peut représenter cette couche comme un graphe de calcul: les entrées et paramètres sont les nœuds sources, les opérations ($\times$, $+$, $\varphi$) sont des transformations, et l’activation $\mathbf{z}$ est le nœud de sortie. Cette perspective sera centrale dans la section sur la dérivation automatique.

Les nœuds en jaune ($W$, $\mathbf{b}$) sont les paramètres (feuilles du graphe); le nœud bleu ($\mathbf{x}$) est l’entrée; le nœud vert ($\mathbf{z}$) est la sortie de la couche.

Rôle de la non-linéarité¶

Sans la fonction d’activation $\varphi$, une couche se réduit à une transformation affine $\mathbf{h} = W\mathbf{x} + \mathbf{b}$. Empiler plusieurs couches linéaires ne fait qu’en produire une autre:

La composition de fonctions linéaires est encore linéaire. Les non-linéarités sont ce qui donne aux réseaux de neurones leur pouvoir expressif.

Fonctions d’activation¶

Nous avons déjà rencontré la sigmoïde en régression logistique au chapitre 3:

Elle transforme un score réel en une valeur dans $(0, 1)$. Sa dérivée est $\sigma'(a) = \sigma(a)(1 - \sigma(a))$, ce qui sera utile pour la rétropropagation. Cependant, la sigmoïde sature pour les grandes valeurs de $|a|$: dans ces régions, la dérivée est proche de zéro.

La tangente hyperbolique est similaire mais centrée autour de zéro:

Ses sorties sont dans $(-1, 1)$. On peut montrer que $\tanh$ est une version recentrée de la sigmoïde. Elle souffre du même problème de saturation.

L’unité linéaire rectifiée (ReLU, de l’anglais rectified linear unit) est aujourd’hui la fonction d’activation la plus utilisée:

Ses avantages sont sa simplicité de calcul et l’absence de saturation pour les valeurs positives. Sa dérivée vaut 1 pour $a > 0$ et 0 pour $a < 0$. Un inconvénient est que les neurones dont la pré-activation est toujours négative ont un gradient nul et cessent d’apprendre: c’est le problème des « neurones morts ».

Plusieurs variantes de ReLU existent pour atténuer ce problème. La Leaky ReLU utilise une petite pente $\alpha \approx 0{,}01$ pour les valeurs négatives: $\text{LeakyReLU}(a) = \max(\alpha a, a)$. La GELU (Gaussian Error Linear Unit), définie par $\text{GELU}(a) = a \cdot \Phi(a)$ où $\Phi$ est la fonction de répartition normale, est utilisée dans les architectures modernes comme les transformeurs.

Source

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
%config InlineBackend.figure_format = 'retina'

a = np.linspace(-4, 4, 400)

sigmoid = lambda x: 1 / (1 + np.exp(-x))
gelu    = lambda x: x * 0.5 * (1 + erf(x / np.sqrt(2)))

d_sigmoid = lambda x: sigmoid(x) * (1 - sigmoid(x))
d_tanh    = lambda x: 1 - np.tanh(x)**2
d_relu    = lambda x: (x > 0).astype(float)
eps = 1e-5
d_gelu    = lambda x: (gelu(x + eps) - gelu(x - eps)) / (2 * eps)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

ax = axes[0]
ax.plot(a, sigmoid(a),       'C0', linewidth=2, label='Sigmoïde')
ax.plot(a, np.tanh(a),       'C1', linewidth=2, label='Tanh')
ax.plot(a, np.maximum(0, a), 'C2', linewidth=2, label='ReLU')
ax.plot(a, gelu(a),          'C3', linewidth=2, label='GELU', linestyle='--')
ax.axhline(0, color='k', linewidth=0.5, linestyle=':')
ax.axvline(0, color='k', linewidth=0.5, linestyle=':')
ax.set_xlabel('$a$ (pré-activation)')
ax.set_ylabel('$\\varphi(a)$')
ax.set_title("Fonctions d'activation")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)

ax = axes[1]
ax.plot(a, d_sigmoid(a), 'C0', linewidth=2, label="$\\sigma'(a)$")
ax.plot(a, d_tanh(a),    'C1', linewidth=2, label="$\\tanh'(a)$")
ax.plot(a, d_relu(a),    'C2', linewidth=2, label="ReLU$'(a)$")
ax.plot(a, d_gelu(a),    'C3', linewidth=2, label="GELU$'(a)$", linestyle='--')
ax.axhline(0, color='k', linewidth=0.5, linestyle=':')
ax.axvline(0, color='k', linewidth=0.5, linestyle=':')
ax.annotate(
    "$\\sigma'(0) = 0{,}25$",
    xy=(0, 0.25), xytext=(1.3, 0.42),
    arrowprops=dict(arrowstyle='->', color='C0', lw=1.5),
    fontsize=9, color='C0'
)
ax.set_xlabel('$a$ (pré-activation)')
ax.set_ylabel("$\\varphi'(a)$")
ax.set_title('Dérivées des fonctions d\'activation')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
ax.set_ylim(-0.1, 1.1)

plt.tight_layout()

La dérivée de la sigmoïde est bornée par 0,25: à chaque couche, le gradient est multiplié par un facteur d’au plus 0,25. Ce phénomène de saturation est la cause principale de la dissolution du gradient dans les réseaux profonds. Nous y reviendrons en détail au chapitre suivant.

Le perceptron multicouche¶

Un perceptron multicouche (MLP, de l’anglais multilayer perceptron) compose plusieurs couches de la forme décrite ci-dessus. Pour un réseau à $L$ couches:

L’entrée $\mathbf{x}$ traverse $L-1$ couches cachées, chacune produisant des activations $\mathbf{z}_\ell$. Ces activations sont les caractéristiques apprises, c’est-à-dire la fonction $\boldsymbol{\phi}(\mathbf{x}; \boldsymbol{\theta}_\phi)$ de l’équation (7). La dernière couche produit la sortie du réseau.

Le graphe de calcul d’un MLP à deux couches cachées est une chaîne d’opérations: chaque couche correspond à une transformation affine suivie d’une non-linéarité.

Ce graphe est une structure de données: il encode toutes les dépendances entre variables. La rétropropagation consiste à le parcourir en sens inverse pour calculer les gradients par rapport à chaque paramètre.

Couche de sortie: le lien avec le maximum de vraisemblance¶

Le traitement de la dernière couche dépend du problème et de notre choix de distribution conditionnelle:

Pour la régression (vraisemblance gaussienne), la couche de sortie est linéaire, sans activation:

La perte est la somme des carrés, $\sum_i (y_i - \mu(\mathbf{x}_i))^2$, cohérente avec l’hypothèse de bruit gaussien.

Pour la classification binaire (vraisemblance de Bernoulli), la couche de sortie applique une sigmoïde:

La perte est l’entropie croisée binaire, exactement comme en régression logistique.

Pour la classification multiclasse (vraisemblance catégorielle), la couche de sortie applique un softmax:

La perte est l’entropie croisée catégorielle.

Un réseau de neurones pour la classification est donc une régression logistique dont les entrées sont des caractéristiques apprises. Les couches cachées construisent une représentation $\mathbf{z}_{L-1} = \boldsymbol{\phi}(\mathbf{x}; \boldsymbol{\theta}_\phi)$ dans laquelle le problème devient (idéalement) linéairement séparable, et la dernière couche effectue la classification linéaire.

Expressivité¶

Un réseau avec une seule couche cachée suffisamment large peut approximer toute fonction continue sur un ensemble compact. Ce résultat, connu sous le nom de théorème d’approximation universelle Hornik et al. (1989), garantit l’expressivité théorique des MLP. Cependant, la largeur requise peut croître exponentiellement avec la complexité de la fonction cible. Les réseaux profonds (avec plusieurs couches) peuvent représenter certaines fonctions de manière beaucoup plus compacte que les réseaux larges mais peu profonds.

La figure ci-dessous illustre cette propriété: un réseau peu profond mais large et un réseau profond mais étroit approximent tous deux la même fonction, mais avec des complexités très différentes.

Source

import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'

np.random.seed(0)

def relu(x):
    return np.maximum(0, x)

def forward_shallow(x, W1, b1, W2, b2):
    """Réseau large: 1 couche cachée, beaucoup de neurones."""
    h = relu(x[:, None] * W1[None, :] + b1[None, :])
    return h @ W2 + b2

def forward_deep(x, params):
    """Réseau profond: plusieurs couches, peu de neurones."""
    h = x[:, None]
    for W, b in params[:-1]:
        h = relu(h @ W + b)
    W, b = params[-1]
    return (h @ W + b).ravel()

# Cible: fonction non triviale
x_grid = np.linspace(0, 1, 200)
f_target = np.sin(2 * np.pi * x_grid) + 0.5 * np.sin(6 * np.pi * x_grid)

# Réseau peu profond, large (1 couche cachée, 40 neurones)
n_wide = 40
W1_s = np.random.randn(n_wide) * 3
b1_s = np.random.randn(n_wide)
# Ajuster W2 par pseudoinverse pour approximer la cible
H_s = relu(x_grid[:, None] * W1_s[None, :] + b1_s[None, :])
W2_s, _, _, _ = np.linalg.lstsq(
    np.column_stack([H_s, np.ones(len(x_grid))]),
    f_target, rcond=None
)
b2_s = W2_s[-1]
W2_s = W2_s[:-1]
pred_shallow = H_s @ W2_s + b2_s

# Réseau profond, étroit (4 couches cachées, 8 neurones)
n_deep = 8
params_deep = []
d_in = 1
for layer in range(4):
    W = np.random.randn(d_in, n_deep) * np.sqrt(2 / d_in)
    b = np.zeros(n_deep)
    params_deep.append((W, b))
    d_in = n_deep
W_out = np.random.randn(d_in, 1) * np.sqrt(2 / d_in)
b_out = np.zeros(1)
params_deep.append((W_out, b_out))

# Ajuster la dernière couche par pseudoinverse
h = x_grid[:, None]
for W, b in params_deep[:-1]:
    h = relu(h @ W + b)
W_out_fit, _, _, _ = np.linalg.lstsq(
    np.column_stack([h, np.ones(len(x_grid))]),
    f_target, rcond=None
)
b_out_fit = W_out_fit[-1]
W_out_fit = W_out_fit[:-1, None]
params_deep[-1] = (W_out_fit, np.array([b_out_fit]))

h = x_grid[:, None]
for W, b in params_deep[:-1]:
    h = relu(h @ W + b)
W_f, b_f = params_deep[-1]
pred_deep = (h @ W_f + b_f).ravel()

# Figure
fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)

for ax, pred, title, n_params in zip(
    axes,
    [pred_shallow, pred_deep],
    [f'Peu profond, large\n(1 couche cachée, {n_wide} neurones)',
     f'Profond, étroit\n(4 couches cachées, {n_deep} neurones)'],
    [n_wide * 2 + n_wide + 1, 1 * n_deep + n_deep + 3 * n_deep * n_deep + n_deep + n_deep + 1]
):
    ax.plot(x_grid, f_target, 'k--', linewidth=2, label='Cible $f(x)$', alpha=0.7)
    ax.plot(x_grid, pred, 'C0', linewidth=2, label='Approximation')
    ax.set_xlabel('$x$')
    ax.set_title(title, fontsize=10)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

axes[0].set_ylabel('$f(x)$')
plt.suptitle("Théorème d'approximation universelle: deux architectures, une même fonction", fontsize=11)
plt.tight_layout()

Les deux architectures approximent raisonnablement bien la même fonction cible. La différence réside dans l’organisation des paramètres: le réseau profond compose des représentations intermédiaires hiérarchiques, ce qui lui permet d’être plus compact pour des fonctions structurées.

À ce stade, vous avez vu la structure d’un réseau de neurones: des couches qui alternent transformations linéaires et non-linéarités, avec une couche de sortie adaptée au problème (régression ou classification). La question suivante est: comment apprendre les paramètres?

Dérivation automatique¶

Pour optimiser les paramètres d’un réseau, nous avons besoin du gradient de la perte par rapport à chaque paramètre. Dans un réseau à $L$ couches, la perte dépend des paramètres de la couche $\ell$ à travers toutes les couches suivantes $\ell+1, \ldots, L$: le calcul du gradient exige d’appliquer la règle de la chaîne à travers tout le graphe de calcul.

Cette section présente la dérivation automatique (DA), également appelée dérivation algorithmique ou, dans la littérature anglophone, automatic differentiation (AD). C’est le cadre général qui formalise ce calcul. Nous commençons par situer la DA parmi les approches de calcul de dérivées, puis nous développons la règle de la chaîne sous forme de produits jacobien-vecteur (JVP et VJP). Nous introduisons ensuite les graphes de calcul (DAG) et formulons les algorithmes du mode avant et du mode arrière sur un graphe arbitraire. La rétropropagation (backpropagation) apparaît alors comme un cas particulier du mode arrière, appliqué au graphe en chaîne d’un réseau de neurones. Enfin, nous montrons comment les bibliothèques modernes (JAX, PyTorch) implémentent ces principes via le traçage d’opérations.

Dérivation numérique, symbolique et automatique¶

Pour calculer la dérivée d’un programme, trois approches existent:

La dérivation numérique approxime la dérivée par différences finies:

Cette méthode est simple à implémenter mais souffre de deux problèmes: elle requiert $O(n)$ évaluations de $f$ pour un gradient en dimension $n$, et elle est sujette aux erreurs d’arrondi (le choix de $\epsilon$ est délicat). Elle reste utile pour vérifier des implémentations de gradient.

La dérivation symbolique applique les règles de dérivation formellement, comme on le ferait à la main. Le résultat est une expression mathématique, pas un nombre. La bibliothèque SymPy permet de s’en convaincre:

import sympy as sp

x = sp.Symbol('x')
f = sp.sin(x**2) * sp.exp(-x)

df = sp.diff(f, x)
print(df)

2*x*exp(-x)*cos(x**2) - exp(-x)*sin(x**2)

L’appel sp.diff(f, x) retourne une nouvelle expression symbolique: $2x\cos(x^2)e^{-x} - \sin(x^2)e^{-x}$. Le système manipule des formules, pas des valeurs numériques. Pour obtenir un nombre, il faut ensuite évaluer cette expression en un point:

df.subs(x, 1.0).evalf()  # évalue la dérivée en x = 1

Cette distinction entre construire une expression et évaluer un nombre est au cœur de la différence entre dérivation symbolique et automatique.

L’approche symbolique produit des résultats exacts, mais les expressions intermédiaires peuvent croître de façon exponentielle. Considérons une composition itérée $h(x) = \sin(\sin(\cdots\sin(x)\cdots))$ sur $k$ niveaux. Chaque application de la règle de la chaîne multiplie l’expression par un facteur $\cos(\cdots)$, et l’expression de la dérivée accumule un produit de cosinus imbriqués dont la taille croît avec $k$. Pour des programmes réels avec des centaines d’opérations, cette croissance rend l’approche impraticable. De plus, la dérivation symbolique requiert que le programme soit représenté sous forme d’expression mathématique, ce qui exclut les structures de contrôle comme les boucles et les conditions.

La dérivation automatique (DA) est une troisième voie. Au lieu de construire l’expression symbolique de la dérivée puis de l’évaluer, elle évalue directement la dérivée en un point donné, en propageant des valeurs numériques à travers le programme. Chaque opération élémentaire (addition, multiplication, $\sin$, $\exp$, ...) est accompagnée de sa règle de dérivation locale, et la règle de la chaîne assemble ces dérivées locales au fur et à mesure de l’exécution. Le résultat est un nombre, la valeur exacte de la dérivée au point considéré, obtenu sans jamais former une expression intermédiaire. Contrairement à la dérivation numérique, ce résultat est exact (aux erreurs de virgule flottante près). Contrairement à la dérivation symbolique, il gère naturellement les boucles et les conditions, puisqu’il opère sur l’exécution concrète du programme.

La rétropropagation n’est rien d’autre que la dérivation automatique en mode arrière, appliquée au programme qui calcule la perte d’un réseau de neurones.

La règle de la chaîne pour les compositions¶

Considérons un réseau comme une composition de fonctions $f = f_L \circ f_{L-1} \circ \cdots \circ f_1$. La jacobienne de cette composition est le produit des jacobiennes individuelles:

où $\mathbf{z}_\ell = f_\ell(\mathbf{z}_{\ell-1})$ sont les valeurs intermédiaires calculées lors de la passe avant. Ce produit de matrices peut être évalué de deux façons, et le choix fait toute la différence.

Produits jacobien-vecteur¶

Le produit $\mathbf{J}_f \cdot \mathbf{v}$ d’une jacobienne par un vecteur peut être calculé sans jamais former la jacobienne complète. Selon la direction de multiplication, on obtient deux opérations distinctes.

Le JVP (Jacobian-Vector Product) propage un vecteur tangent $\mathbf{v}$ de gauche à droite:

Chaque étape multiplie une jacobienne locale par un vecteur, ce qui coûte $O(mn)$ au lieu de $O(m n^2)$ pour le produit par une matrice. Le calcul se fait dans le même sens que la passe avant: c’est le mode avant de la dérivation automatique.

Le VJP (Vector-Jacobian Product) propage un vecteur adjoint $\mathbf{u}^\top$ de droite à gauche:

Le calcul se fait dans le sens inverse de la passe avant: c’est le mode arrière.

Pour une perte scalaire $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$, le gradient $\nabla_\mathbf{x} \mathcal{L}$ est exactement un VJP avec $\mathbf{u} = 1$. Le mode arrière calcule donc le gradient par rapport à tous les paramètres en une seule passe arrière, quel que soit le nombre de paramètres. C’est pourquoi la rétropropagation utilise le mode arrière.

La figure ci-dessous contraste les deux modes sur une chaîne de trois fonctions $f_1 \circ f_2 \circ f_3$. Le mode avant (JVP) propage un vecteur tangent de gauche à droite, ce qui coûte une passe par paramètre. Le mode arrière (VJP) propage l’adjoint de droite à gauche en une seule passe.

Mode avant (JVP): le vecteur tangent $\tilde{v}$ se propage de gauche à droite.

Mode arrière (VJP): l’adjoint $\bar{u}$ se propage de droite à gauche.

Pour une perte scalaire avec $n$ paramètres, le mode avant nécessite $n$ passes (une par direction de base $\mathbf{e}_i$), tandis que le mode arrière calcule tout le gradient en une seule passe. C’est l’argument central qui justifie la rétropropagation dans les réseaux avec des millions de paramètres.

Le carnet Produits jacobien-vecteur en mode inverse (VJP) illustre ce mécanisme pas à pas sur un réseau à trois couches, en vérifiant à chaque étape que les VJP produisent le même résultat que les produits matriciels explicites.

Graphes de calcul¶

Les deux modes ont été présentés pour une composition en chaîne. Mais les programmes réels ont des structures plus riches: une variable peut intervenir dans plusieurs opérations, créant des embranchements dans le graphe de calcul.

Toute expression arithmétique peut se décomposer en une séquence d’opérations élémentaires. Pour $f(x, y) = \sin(x) \cdot (x + y)$, cette décomposition introduit deux variables intermédiaires:

On représente ces dépendances par un graphe orienté acyclique (DAG): chaque noeud est une valeur (entrée, intermédiaire ou sortie) et chaque arête indique qu’une valeur est utilisée pour calculer une autre.

Règle de la chaîne sur un DAG¶

Remarquez que $x$ a deux arêtes sortantes: il alimente $v_1 = \sin(x)$ et $v_2 = x + y$. Cela signifie que $x$ contribue à $f$ par deux chemins distincts dans le graphe. Notons $\phi_1(x) = \sin(x)$, $\phi_2(x,y) = x + y$ et $\phi_3(v_1,v_2) = v_1 \cdot v_2$ les opérations locales de chaque noeud. La fonction composée est $f(x,y) = \phi_3(\phi_1(x),\, \phi_2(x,y))$. En appliquant la règle de la chaîne multivariée, la dérivée totale de $f$ par rapport à $x$ est:

La somme comporte deux termes, un par chemin de $x$ à $f$ dans le graphe:

• Chemin $x \to v_1 \to f$: le produit des jacobiennes locales le long des arêtes, $D_{v_1}\phi_3 \cdot D_x \phi_1$

• Chemin $x \to v_2 \to f$: le produit $D_{v_2}\phi_3 \cdot D_x \phi_2$

C’est la structure générale: la dérivée totale par rapport à une variable est la somme sur tous les chemins de cette variable à la sortie, où chaque chemin contribue le produit des jacobiennes locales le long de ses arêtes. Un noeud avec $k$ arêtes sortantes génère $k$ termes dans cette somme.

De manière générale, notons $\text{pred}(v)$ l’ensemble des prédécesseurs d’un noeud $v$ (les noeuds dont $v$ dépend directement) et $\text{succ}(u)$ l’ensemble de ses successeurs (les noeuds qui dépendent directement de $u$). La règle de la chaîne s’écrit dans les deux sens:

Direction avant (tangentes). Étant donné un vecteur tangent $\dot{x}_i$ pour chaque entrée, la tangente d’un noeud intermédiaire se propage vers l’avant. Le noeud $v$ reçoit les tangentes de tous ses prédécesseurs et les combine:

Chaque arête entrante contribue un JVP (jacobienne locale $\times$ tangente du prédécesseur). Le noeud $v$ somme ces contributions.

Direction arrière (adjoints). Étant donné un adjoint $\bar{f} = 1$ pour la sortie, l’adjoint de chaque noeud se propage vers l’arrière. Le noeud $u$ reçoit les adjoints de tous ses successeurs:

Chaque arête sortante (parcourue à rebours) contribue un VJP (adjoint du successeur $\times$ jacobienne locale). Le noeud $u$ somme ces contributions.

Les deux formules sont symétriques: la première applique $D_u \phi_v$ à droite d’un vecteur tangent (JVP), la seconde applique $D_u \phi_v$ à gauche d’un vecteur adjoint (VJP). C’est la règle de la chaîne multivariée.

Tri topologique¶

Les deux formules ci-dessus posent un problème d’ordre. Pour calculer $\dot{v}$ (direction avant), il faut d’abord connaître $\dot{u}$ pour tous les prédécesseurs $u$ de $v$. Pour calculer $\bar{u}$ (direction arrière), il faut d’abord connaître $\bar{v}$ pour tous les successeurs $v$ de $u$.

Un tri topologique fournit un ordre de traitement qui respecte ces contraintes: chaque noeud apparaît après tous ses prédécesseurs. La passe avant suit cet ordre; la passe arrière le parcourt à rebours.

Le diagramme ci-dessous montre un tri topologique valide pour notre exemple. Les numéros indiquent l’ordre de traitement de la passe arrière (qui parcourt la liste à rebours):

Remarquez que $x$ est traité en dernier (⑥) par la passe arrière: comme $x$ contribue à deux branches ($v_1$ et $v_2$), il faut avoir accumulé les deux contributions avant de pouvoir calculer $\bar{x}$.

L’algorithme classique de tri topologique repose sur un parcours en profondeur (DFS):

Mode avant (forward-mode AD)¶

Le mode avant propage les tangentes dans le même sens que l’exécution du programme, de l’entrée vers la sortie. On calcule simultanément la valeur de chaque noeud et sa tangente.

Une passe avant calcule un seul produit jacobien-vecteur $\mathbf{J}_f \dot{\mathbf{x}}$, c’est-à-dire la dérivée directionnelle de $f$ dans la direction $\dot{\mathbf{x}}$. Pour obtenir le gradient complet d’une fonction scalaire $f: \mathbb{R}^n \to \mathbb{R}$, il faudrait effectuer $n$ passes (une par vecteur de base $\mathbf{e}_i$). Le mode avant est donc efficace quand le nombre d’entrées est petit par rapport au nombre de sorties.

Mode arrière (reverse-mode AD)¶

Le mode arrière procède en deux temps: d’abord une passe avant qui calcule et stocke les valeurs de chaque noeud, puis une passe arrière qui propage les adjoints de la sortie vers les entrées.

Une seule passe arrière calcule le gradient par rapport à toutes les entrées. Pour une perte scalaire avec $n$ paramètres, le mode arrière calcule le gradient complet en une passe, là où le mode avant en nécessiterait $n$. C’est l’argument central qui justifie l’usage du mode arrière pour l’entraînement des réseaux de neurones.

Le mode arrière a été décrit pour la première fois par Seppo Linnainmaa Linnainmaa (1970) dans sa thèse de maîtrise à l’Université d’Helsinki. Dans le contexte de l’apprentissage profond, cet algorithme porte le nom de rétropropagation (backpropagation), popularisé par Rumelhart, Hinton et Williams en 1986 Rumelhart et al. (1986). Un réseau à $K$ couches définit un DAG en chaîne $f_1 \circ f_2 \circ \cdots \circ f_K$, où chaque noeud n’a qu’un seul prédécesseur et un seul successeur. La passe arrière se simplifie alors en une boucle sur les couches $K, K-1, \ldots, 1$. Mais le mode arrière est plus général: il s’applique à tout programme différentiable, y compris ceux qui contiennent des embranchements, des variables réutilisées, des boucles ou des conditions. C’est ce qui permet à des bibliothèques comme JAX ou PyTorch de différentier n’importe quelle fonction Python.

L’animation interactive ci-dessous illustre les deux algorithmes sur le DAG de $f(x,y) = \sin(x) \cdot (x + y)$ avec les valeurs $(x, y) = (0.5, \; 1.2)$. En mode avant, les tangentes se propagent de gauche à droite; en mode arrière, les adjoints se propagent de droite à gauche. Observez en particulier comment $\bar{x}$ accumule les contributions de ses deux successeurs $\phi_1$ et $\phi_2$ lors de la passe arrière.

from IPython.display import HTML
from pathlib import Path
import html as _html

_content = Path("_static/ad_stepper.html").read_text()
_doc = f'<!DOCTYPE html><html><head><meta charset="utf-8"></head><body style="margin:0;padding:0">{_content}</body></html>'
HTML(f'<iframe srcdoc="{_html.escape(_doc, quote=True)}" width="100%" height="580" style="border:none;" scrolling="no"></iframe>')

/Users/pierre-luc.bacon/Documents/mlbook/.venv/lib/python3.12/site-packages/IPython/core/display.py:447: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")

Règles VJP: une bibliothèque d’opérateurs adjoints¶

L’algorithme du mode arrière fait intervenir les jacobiennes locales $D_u \phi_v$ à chaque étape. Une question reste ouverte: comment calcule-t-on $\bar{v} \, D_u \phi_v$ efficacement, sans construire la jacobienne complète?

La réponse repose sur une distinction fondamentale. La jacobienne $\mathbf{J}_f(\mathbf{x}) \in \mathbb{R}^{m \times n}$ est une représentation matricielle d’un objet plus abstrait: la différentielle $df_\mathbf{x}$, qui est un opérateur linéaire $df_\mathbf{x}: \mathbb{R}^n \to \mathbb{R}^m$. Ce qui importe dans le mode arrière n’est pas $df_\mathbf{x}$ lui-même, mais son opérateur adjoint $df_\mathbf{x}^*: \mathbb{R}^m \to \mathbb{R}^n$, qui envoie les vecteurs du co-domaine vers le domaine (c’est-à-dire qui propage le signal en sens inverse). En coordonnées, le VJP est $\mathbf{u}^\top \mathbf{J}_f(\mathbf{x})$: le produit d’un vecteur adjoint (ligne) par la jacobienne.

Pour définir un opérateur linéaire, il n’est pas nécessaire d’en donner la matrice: on peut spécifier son action sur des vecteurs. Une bibliothèque de DA (JAX, PyTorch) maintient pour chaque opération primitive une règle VJP: une fonction qui calcule directement $\mathbf{u}^\top \mathbf{J}_f(\mathbf{x})$ à partir de $\mathbf{u}$, $\mathbf{x}$, et éventuellement $f(\mathbf{x})$, en n’utilisant que des opérations arithmétiques simples.

Lorsque deux opérations se composent, $h = g \circ f$, la règle VJP de $h$ est le produit des règles VJP de $g$ et $f$:

La passe arrière n’est rien d’autre que l’exécution récursive de ces règles, de la sortie vers l’entrée. Le système est entièrement sans matrice jacobienne explicite: aucune matrice $\mathbf{J}_f$ n’est jamais construite ni stockée.

Le tableau ci-dessous liste les règles VJP pour les opérations clés d’un MLP. À chaque fois, la règle VJP évite de former la jacobienne correspondante:

L’exemple de l’activation élémentaire illustre parfaitement le bénéfice: la jacobienne serait une matrice $m \times m$ coûtant $O(m^2)$ en mémoire, alors que la règle VJP, $\mathbf{u} \odot \varphi'(\mathbf{a})$, est un produit élément par élément en $O(m)$.

En JAX, jax.custom_vjp permet d’enregistrer exactement ce type de règle pour une opération personnalisée. Écrire une règle VJP correcte pour un nouvel opérateur est une compétence essentielle en apprentissage profond avancé; les exercices 9 à 11 vous entraînent à cette dérivation.

Exemple: MLP avec une couche cachée¶

Prenons un réseau à une couche cachée avec la perte des moindres carrés:

Le graphe de calcul de ce réseau rend explicites toutes les dépendances. Les nœuds en jaune sont les paramètres (feuilles du graphe); la passe avant suit les flèches de gauche à droite, et la passe arrière les remonte de droite à gauche.

La passe arrière calcule les adjoints en remontant ce graphe nœud par nœud, en appliquant les règles VJP de chaque opération.

La passe avant calcule les valeurs intermédiaires:

La passe arrière propage les adjoints en sens inverse, couche par couche. La notation $\bar{v}$ désigne l’adjoint du noeud $v$, c’est-à-dire la sensibilité de $\mathcal{L}$ à $v$:

où $\odot$ désigne le produit élément par élément. Chaque ligne utilise uniquement des quantités déjà calculées, soit lors de la passe avant ($\mathbf{z}_1$, $\mathbf{a}_1$, $\mathbf{x}$), soit lors des étapes précédentes de la passe arrière. La structure est toujours la même: l’adjoint des pré-activations d’une couche est propagé vers l’arrière pour obtenir l’adjoint de la couche précédente.

Point de contrôle: Si vous pouvez suivre cet exemple du début à la fin, vous avez compris le mécanisme de la rétropropagation. C’est une instance de l’algorithme du mode arrière de la section précédente, spécialisée à un réseau en chaîne. Si certaines étapes restent floues, l’exercice 3 vous permettra de refaire ce calcul vous-même avec des valeurs numériques.

L’animation interactive ci-dessous déroule la passe avant et la rétropropagation sur un réseau à deux couches avec des valeurs numériques concrètes. Observez comment chaque couche produit les gradients de ses paramètres ($\bar{w}$, $\bar{b}$) tout en propageant le signal vers l’arrière, et comment la dérivée de l’activation ($\sigma'$) joue le rôle de «porte» qui laisse passer ou bloque le gradient.

from IPython.display import HTML
from pathlib import Path
import html as _html

_content = Path("_static/nn_backprop_stepper.html").read_text()
_doc = f'<!DOCTYPE html><html><head><meta charset="utf-8"></head><body style="margin:0;padding:0">{_content}</body></html>'
HTML(f'<iframe srcdoc="{_html.escape(_doc, quote=True)}" width="100%" height="580" style="border:none;" scrolling="no"></iframe>')

La liste de Wengert¶

En 1964, Wengert Wengert (1964) a proposé de représenter toute fonction calculable comme une séquence ordonnée d’opérations élémentaires, chacune ayant une dérivée connue. Cette séquence, appelée liste de Wengert (ou tape, bande), est le graphe de calcul sérialisé en ordre topologique. La contribution de Wengert était de rendre la dérivation algorithmique: en décomposant un programme en étapes atomiques et en les écrivant dans l’ordre, une machine peut appliquer la règle de la chaîne mécaniquement, sans intervention humaine.

En mode avant, cette liste guide la propagation des tangentes: chaque étape calcule simultanément sa valeur et sa dérivée, sans rien retenir. Mais en mode arrière, la bande devient indispensable comme structure de stockage: il faut rejouer les opérations à rebours, ce qui exige de conserver les valeurs intermédiaires de la passe avant. C’est ce rôle de stockage qui domine dans les bibliothèques modernes de DA (JAX, PyTorch), puisque l’entraînement des réseaux utilise le mode arrière.

Concrètement, pendant la passe avant, chaque opération enregistre sur la bande ses entrées, sa sortie, et une fonction VJP locale. Cette fonction est construite au moment de l’opération, et elle capture les valeurs intermédiaires dont elle aura besoin plus tard pour calculer le gradient. En programmation, on appelle cela une fermeture (closure). Prenons l’opération $v_3 = v_1 \cdot v_2$ comme exemple. Au moment du calcul, la passe avant connaît les valeurs de $v_1$ et $v_2$. La règle VJP de la multiplication est $(\bar{v}_1, \bar{v}_2) = (v_2 \cdot \bar{v}_3,\; v_1 \cdot \bar{v}_3)$: elle a besoin des valeurs $v_1$ et $v_2$ de la passe avant. La fermeture les capture:

# Pendant la passe avant, au moment de v3 = v1 * v2 :
v1_val, v2_val = v1, v2           # valeurs connues maintenant

def mul_vjp(v3_bar):              # sera appelée pendant la passe arrière
    return (v2_val * v3_bar,      # gradient pour v1
            v1_val * v3_bar)      # gradient pour v2

tape.append(mul_vjp)              # on stocke la fermeture, pas les valeurs brutes

La fonction mul_vjp ne sera appelée que plus tard, pendant la passe arrière, avec l’adjoint $\bar{v}_3$ comme argument. Mais elle a déjà accès à v1_val et v2_val, capturés au moment de sa création. C’est ce mécanisme de capture qui rend la bande autonome: chaque entrée contient tout ce qu’il faut pour calculer sa contribution au gradient, sans consulter à nouveau le programme original.

La passe arrière rejoue la liste à rebours, en appelant chaque VJP locale dans l’ordre inverse. Contrairement au DAG, la bande est une liste ordonnée (un tableau linéaire), ce qui la rend simple à parcourir dans les deux sens. Les flèches pleines montrent l’enregistrement (passe avant); les flèches pointillées montrent le rejeu (passe arrière).

La bande démarre vide et grandit à chaque opération (flèches pleines, gauche à droite). Une fois $f$ calculée, la passe arrière initialise $\bar{v}_3 = 1$ puis remonte la liste à rebours (flèches pointillées, droite à gauche).

Le tableau ci-dessous détaille le contenu de chaque entrée et les formules de VJP associées. La barre $\bar{v}$ désigne l’adjoint $\partial f / \partial v$; la passe arrière parcourt les étapes 3 → 2 → 1.

Chaque ligne de la bande correspond à une opération élémentaire. La passe arrière part de l’adjoint $\bar{v}_3 = 1$ (gradient de $f$ par rapport à lui-même) et remonte: l’étape 3 envoie des gradients à $v_1$ et $v_2$, puis les étapes 2 et 1 envoient leurs gradients à $x$ et $y$.

Le traceur¶

Comment une bibliothèque comme JAX construit-elle cette bande automatiquement, sans modifier le programme utilisateur? La réponse est le traceur (tracer).

Lorsque JAX différentie une fonction, il ne l’appelle pas avec des nombres ordinaires. Il l’appelle avec des objets spéciaux, des traceurs, qui se font passer pour des nombres mais enregistrent discrètement toutes les opérations qu’on leur applique.

Concrètement, un traceur est un objet qui:

1. Stocke une valeur concrète (le résultat numérique de l’opération),

2. Enregistre l’opération sur la bande, avec ses entrées et une fermeture (closure) qui sait calculer les gradients locaux,

3. Retourne un nouveau traceur comme résultat, de sorte que les opérations suivantes soient également interceptées.

Le diagramme ci-dessous montre comment les objets traceurs se construisent et se connectent lors de l’évaluation de $f(x,y) = \sin(x) \cdot (x + y)$. Chaque noeud porte sa valeur concrète et son adjoint (initialement nul); les arêtes représentent les dépendances enregistrées par chaque fermeture.

La passe arrière initialise $\bar{v}_3 = 1$, puis parcourt les arêtes à rebours: chaque fermeture accumule les adjoints dans les noeuds parents. Quand les deux fermetures de $v_1$ et $v_2$ ont été appelées, x.grad contient la somme des deux contributions.

Le tableau ci-dessous montre la correspondance entre l’exécution Python et la bande construite automatiquement. Chaque ligne de code qui effectue une opération tracée ajoute une entrée sur la bande, avec la règle VJP correspondante.

La passe arrière parcourt la bande de bas en haut (étapes 3 → 2 → 1), en appelant chaque règle VJP.

L’exécution Python se déroule normalement, ligne par ligne. Python ne sait pas qu’il trace un graphe: il appelle simplement les méthodes __add__, __mul__, sin sur les objets traceurs, et ces méthodes enregistrent discrètement les opérations. Quand l’exécution est terminée, la bande est complète, et la passe arrière peut s’exécuter.

Ce mécanisme explique pourquoi la dérivation automatique gère naturellement les boucles et les conditions: Python les exécute normalement, et les traceurs enregistrent les opérations qui sont effectivement effectuées lors de cette exécution particulière.

L’astuce d’importation¶

Le mécanisme du traceur explique une convention qui surprend souvent les débutants. Dans tout code JAX, on écrit:

import jax.numpy as jnp   # et non: import numpy as np

Pourquoi? Lorsque JAX différentie une fonction, il lui passe des traceurs à la place des tableaux NumPy ordinaires. Si vous appelez np.sin(tracer), NumPy ne connaît pas les traceurs: il va tenter de convertir l’objet en tableau numérique, ce qui casse la trace et donne un résultat incorrect (ou lève une erreur).

En revanche, jnp.sin(tracer) est une opération que JAX connaît. JAX intercepte l’appel, enregistre l’opération sur la bande, calcule la valeur concrète, et retourne un nouveau traceur. La trace reste intacte.

import jax
import jax.numpy as jnp
import numpy as np

def f_jnp(x):
    return jnp.sin(x) * x  # correct: jnp intercepte le traceur

def f_np(x):
    return np.sin(x) * x   # incorrect: np ne comprend pas les traceurs

grad_jnp = jax.grad(f_jnp)(1.0)   # fonctionne: retourne cos(1)*1 + sin(1) ≈ 1.382
# grad_np = jax.grad(f_np)(1.0)   # lèverait une erreur ou donnerait un résultat faux

jax.numpy est un espace de noms qui réimplémente toutes les fonctions NumPy de manière à intercepter les traceurs. Pour les tableaux NumPy ordinaires (sans traceur), jnp et np produisent les mêmes résultats numériques. La différence n’apparaît que pendant la trace.

Implémentation minimale¶

Cette sous-section est optionnelle pour IFT3395. Elle montre comment implémenter un moteur de dérivation automatique en mode arrière en une soixantaine de lignes de Python pur.

La section précédente a identifié trois mécanismes: les règles VJP locales, la bande d’enregistrement, et la passe arrière. Nous allons maintenant les assembler en une implémentation fonctionnelle, structurée comme le seraient JAX ou autograd en version simplifiée Maclaurin et al. (2015). L’architecture se décompose en trois parties:

1. Une bibliothèque de règles VJP. Pour chaque opération primitive, une fonction qui prend les résidus (valeurs de la passe avant nécessaires au calcul du gradient) et le cotangent amont $\bar{v}$, et retourne les cotangents pour chaque entrée. Ces fonctions sont les mêmes que celles du tableau de la section précédente.

2. Un traceur avec bande. Un objet Var qui encapsule un flottant et enregistre chaque opération sur une bande globale. C’est l’analogue simplifié des traceurs de JAX.

3. Une fonction grad. Un opérateur d’ordre supérieur, analogue à jax.grad, qui trace la fonction puis parcourt la bande à rebours en appelant les règles VJP.

import math

# ---- 1. Bibliothèque de règles VJP ----
# Signature commune: vjp(résidus, cotangent_sortie) → cotangents_entrées

def add_vjp(res, g):
    return (g, g)                        # ∂(a+b)/∂a = 1, ∂(a+b)/∂b = 1

def mul_vjp(res, g):
    a, b = res
    return (b * g, a * g)                # ∂(a·b)/∂a = b, ∂(a·b)/∂b = a

def sin_vjp(res, g):
    (a,) = res
    return (math.cos(a) * g,)            # ∂sin(a)/∂a = cos(a)

def relu_vjp(res, g):
    (a,) = res
    return (float(a > 0) * g,)           # ∂relu(a)/∂a = 𝟙(a > 0)


# ---- 2. Traceur et bande ----

_tape = []      # bande globale: [(vjp_fn, résidus, ids_entrées, id_sortie)]
_n_vars = 0     # compteur d'identifiants

class Var:
    """Traceur scalaire: encapsule un flottant et un identifiant unique."""

    def __init__(self, data):
        global _n_vars
        self.data = float(data)
        self.id = _n_vars
        _n_vars += 1

    def _record(self, vjp_fn, res, inputs, out_data):
        """Enregistre une opération sur la bande et retourne un nouveau Var."""
        out = Var(out_data)
        _tape.append((vjp_fn, res, [v.id for v in inputs], out.id))
        return out

    def __add__(self, other):
        other = other if isinstance(other, Var) else Var(other)
        return self._record(add_vjp, (self.data, other.data),
                            [self, other], self.data + other.data)

    def __radd__(self, other): return self.__add__(other)

    def __mul__(self, other):
        other = other if isinstance(other, Var) else Var(other)
        return self._record(mul_vjp, (self.data, other.data),
                            [self, other], self.data * other.data)

    def __rmul__(self, other): return self.__mul__(other)

    def sin(self):
        return self._record(sin_vjp, (self.data,),
                            [self], math.sin(self.data))

    def relu(self):
        return self._record(relu_vjp, (self.data,),
                            [self], max(0.0, self.data))


# ---- 3. Fonction grad (analogue à jax.grad) ----

def grad(f):
    """Retourne une fonction qui calcule le gradient de f."""
    def grad_fn(*args):
        global _tape, _n_vars
        _tape, _n_vars = [], 0               # réinitialiser la bande

        # Passe avant: tracer l'exécution
        traced = [Var(a) for a in args]
        result = f(*traced)

        # Passe arrière: propager les cotangents
        adjoints = [0.0] * _n_vars
        adjoints[result.id] = 1.0            # ∂f/∂f = 1

        for vjp_fn, res, in_ids, out_id in reversed(_tape):
            cotangents = vjp_fn(res, adjoints[out_id])
            for idx, ct in zip(in_ids, cotangents):
                adjoints[idx] += ct          # accumulation (embranchement)

        return tuple(adjoints[v.id] for v in traced)
    return grad_fn

La séparation en trois parties n’est pas un choix esthétique: c’est la structure réelle des bibliothèques de DA. Dans JAX, les règles VJP sont enregistrées via jax.custom_vjp, la bande est construite par le traceur interne, et jax.grad orchestre la passe arrière. Notre implémentation reproduit cette architecture en miniature.

Vérifions sur $f(x, y) = \sin(x) \cdot (x + y)$:

import math

# --- Valeurs de test ---
x0, y0 = 0.5, 1.2

# --- Avec notre moteur de DA ---
def f(x, y):
    return x.sin() * (x + y)

df_dx, df_dy = grad(f)(x0, y0)

print(f'f({x0}, {y0})          = {math.sin(x0) * (x0 + y0):.6f}')
print(f'∂f/∂x (AD)           = {df_dx:.6f}')
print(f'∂f/∂y (AD)           = {df_dy:.6f}')

# --- Vérification analytique ---
df_dx_exact = math.cos(x0) * (x0 + y0) + math.sin(x0)
df_dy_exact = math.sin(x0)
print(f'∂f/∂x (exact)        = {df_dx_exact:.6f}')
print(f'∂f/∂y (exact)        = {df_dy_exact:.6f}')

# --- Vérification par différences finies ---
eps = 1e-5
df_dx_num = (math.sin(x0+eps)*(x0+eps+y0) - math.sin(x0-eps)*(x0-eps+y0)) / (2*eps)
df_dy_num = (math.sin(x0)*(x0+y0+eps)     - math.sin(x0)*(x0+y0-eps))     / (2*eps)
print(f'∂f/∂x (diff. fin.)  = {df_dx_num:.6f}')
print(f'∂f/∂y (diff. fin.)  = {df_dy_num:.6f}')

f(0.5, 1.2)          = 0.815023
∂f/∂x (AD)           = 1.971316
∂f/∂y (AD)           = 0.479426
∂f/∂x (exact)        = 1.971316
∂f/∂y (exact)        = 0.479426
∂f/∂x (diff. fin.)  = 1.971316
∂f/∂y (diff. fin.)  = 0.479426

Les trois méthodes sont en accord. Remarquez que grad est une fonction d’ordre supérieur qui retourne une nouvelle fonction, exactement comme jax.grad. L’accumulation des adjoints (ligne adjoints[idx] += ct) gère automatiquement le cas où une variable contribue à plusieurs branches du calcul: c’est la somme des deux chemins pour $x$.

La programmation différentiable¶

Les bibliothèques modernes comme JAX, PyTorch et TensorFlow implémentent la dérivation automatique de manière générale: toute fonction composée d’opérations dont on connaît les dérivées locales peut être différentiée automatiquement. C’est le paradigme de la programmation différentiable (differentiable programming).

Au lieu de dériver manuellement les gradients pour chaque architecture, nous écrivons la passe avant comme un programme ordinaire, et la bibliothèque se charge de calculer les gradients.

Voici un exemple avec JAX. Nous définissons la passe avant d’un MLP à une couche cachée, puis utilisons jax.grad pour obtenir automatiquement la fonction qui calcule les gradients:

import jax
import jax.numpy as jnp

def predict(params, x):
    """Passe avant d'un MLP à une couche cachée."""
    W1, b1, W2, b2 = params
    h = jnp.tanh(W1 @ x + b1)  # couche cachée
    return W2 @ h + b2           # couche de sortie

def loss_fn(params, x, y):
    """Perte des moindres carrés."""
    y_pred = predict(params, x)
    return 0.5 * jnp.sum((y_pred - y) ** 2)

# jax.grad retourne une FONCTION qui calcule le gradient
grad_fn = jax.grad(loss_fn)

# Un seul appel donne les gradients par rapport à tous les paramètres
grads = grad_fn(params, x, y)

La fonction loss_fn est un programme Python ordinaire. L’appel jax.grad(loss_fn) produit une nouvelle fonction qui calcule le gradient par rapport au premier argument (params). Aucune dérivation manuelle n’est nécessaire: JAX applique la règle de la chaîne automatiquement, en mode arrière, sur la trace d’exécution du programme.

Ce paradigme change la façon de penser les modèles. Au lieu de concevoir une architecture puis de dériver ses gradients, on conçoit un programme de calcul quelconque, avec des boucles, des conditions, des appels de fonctions, et on le différentie automatiquement. La seule contrainte est que les opérations soient différentiables (ou différentiables presque partout, comme ReLU).

Implémentation¶

Cette section réunit les concepts du chapitre dans une implémentation complète d’un MLP en NumPy. L’optimiseur Adam utilisé ici est décrit en détail au chapitre suivant. Le code est volontairement auto-contenu et commenté pas à pas: l’objectif est de rendre le lien entre les équations et le code aussi direct que possible.

Classe MLP avec Adam¶

import numpy as np

class MLP:
    """
    Perceptron multicouche à une couche cachée.
    Activation cachée: ReLU. Activation de sortie: sigmoïde (classification binaire).
    Optimiseur: Adam.
    """

    def __init__(self, n_input, n_hidden, n_output,
                 eta=1e-3, beta1=0.9, beta2=0.999, eps=1e-8,
                 lam=0.0, p_drop=0.0, seed=0):
        rng = np.random.default_rng(seed)
        # Initialisation He pour ReLU
        self.W1 = rng.standard_normal((n_input,  n_hidden)) * np.sqrt(2 / n_input)
        self.b1 = np.zeros(n_hidden)
        # Initialisation Glorot pour la couche de sortie (sigmoïde)
        self.W2 = rng.standard_normal((n_hidden, n_output)) * np.sqrt(2 / (n_hidden + n_output))
        self.b2 = np.zeros(n_output)

        self.eta   = eta
        self.beta1 = beta1;  self.beta2 = beta2;  self.eps = eps
        self.lam   = lam     # décroissance des poids
        self.p_drop = p_drop # taux de dropout

        # État interne Adam (moments)
        self._t = 0
        self._m = {k: np.zeros_like(v)
                   for k, v in [('W1',self.W1),('b1',self.b1),
                                 ('W2',self.W2),('b2',self.b2)]}
        self._s = {k: np.zeros_like(v) for k, v in self._m.items()}

    # ------------------------------------------------------------------
    def _relu(self, x):     return np.maximum(0, x)
    def _sigmoid(self, x):  return 1 / (1 + np.exp(-np.clip(x, -50, 50)))

    # ------------------------------------------------------------------
    def forward(self, X, training=False):
        """Passe avant. Retourne les sorties et les caches."""
        a1 = X @ self.W1 + self.b1          # pré-activations couche cachée
        z1 = self._relu(a1)                  # activations ReLU

        # Dropout (entraînement uniquement)
        if training and self.p_drop > 0:
            mask = (np.random.rand(*z1.shape) > self.p_drop) / (1 - self.p_drop)
            z1 = z1 * mask
        else:
            mask = np.ones_like(z1)

        a2   = z1 @ self.W2 + self.b2       # pré-activations couche de sortie
        pred = self._sigmoid(a2)             # probabilités
        cache = {'X': X, 'a1': a1, 'z1': z1, 'mask': mask}
        return pred, cache

    # ------------------------------------------------------------------
    def backward(self, pred, y, cache):
        """
        Passe arrière. Retourne les gradients par rapport à tous les paramètres.
        y: vecteur colonne de cibles binaires.
        """
        B = len(y)
        X, a1, z1, mask = cache['X'], cache['a1'], cache['z1'], cache['mask']

        # Gradient de l'entropie croisée + sigmoïde: dp = pred - y
        dp  = (pred - y) / B

        # Couche de sortie
        dW2 = z1.T @ dp + self.lam * self.W2 / B
        db2 = dp.sum(axis=0)

        # Propagation vers la couche cachée
        dz1 = dp @ self.W2.T
        dz1 = dz1 * mask            # rétropropagation à travers le dropout
        da1 = dz1 * (a1 > 0)        # dérivée de ReLU

        # Couche cachée
        dW1 = X.T @ da1 + self.lam * self.W1 / B
        db1 = da1.sum(axis=0)

        return {'W1': dW1, 'b1': db1, 'W2': dW2, 'b2': db2}

    # ------------------------------------------------------------------
    def _adam_update(self, grads):
        """Applique une étape Adam à tous les paramètres."""
        self._t += 1
        for name, param in [('W1',self.W1),('b1',self.b1),
                              ('W2',self.W2),('b2',self.b2)]:
            g = grads[name]
            self._m[name] = self.beta1 * self._m[name] + (1 - self.beta1) * g
            self._s[name] = self.beta2 * self._s[name] + (1 - self.beta2) * g**2
            mhat = self._m[name] / (1 - self.beta1**self._t)
            shat = self._s[name] / (1 - self.beta2**self._t)
            param -= self.eta * mhat / (np.sqrt(shat) + self.eps)

    # ------------------------------------------------------------------
    def train_step(self, X, y):
        """Une étape d'entraînement sur un mini-lot (X, y)."""
        pred, cache = self.forward(X, training=True)
        grads = self.backward(pred, y.reshape(-1,1).astype(float), cache)
        self._adam_update(grads)
        y_col = y.reshape(-1,1).astype(float)
        p = np.clip(pred, 1e-7, 1-1e-7)
        loss = -np.mean(y_col*np.log(p) + (1-y_col)*np.log(1-p))
        return loss

    def predict_proba(self, X):
        pred, _ = self.forward(X, training=False)
        return pred

    def predict(self, X):
        return (self.predict_proba(X) >= 0.5).astype(int).ravel()

Le MLP en pratique¶

Les sections précédentes ont présenté le MLP comme un objet mathématique: une composition de transformations affines et de non-linéarités. Mais à quoi ressemble cette composition pour un problème concret de régression ou de classification? La réponse nous ramène directement aux modèles linéaires des chapitres 2 et 3.

Régression avec un MLP¶

Au chapitre 2, la régression linéaire prédit la moyenne d’une gaussienne conditionnelle par une transformation affine:

Passer à un MLP revient à composer des transformations affines et des non-linéarités avant cette sortie linéaire. Pour un réseau à deux couches cachées avec $\mathbf{x} \in \mathbb{R}^d$:

où $W_1 \in \mathbb{R}^{h \times d}$, $W_2 \in \mathbb{R}^{h \times h}$, $\mathbf{w}_3 \in \mathbb{R}^h$, et $\varphi$ est une activation (typiquement ReLU). La dernière couche est linéaire (pas d’activation), et la perte reste la somme des carrés, exactement comme au chapitre 2.

Si l’on veut prédire un vecteur $\mathbf{y} \in \mathbb{R}^K$ (par exemple des coordonnées, ou plusieurs cibles simultanément), la dernière couche devient une transformation affine $W_3 \in \mathbb{R}^{K \times h}$:

La perte est la somme des carrés sur les $K$ sorties: $\sum_{k=1}^K (y_k - \hat{y}_k)^2$.

Classification avec un MLP¶

Au chapitre 3, la régression logistique modélise la probabilité de la classe positive par:

Un MLP pour la classification binaire compose des couches cachées avant cette sigmoïde:

La perte est l’entropie croisée binaire, comme en régression logistique. Pour la classification multiclasse avec $K$ classes, la sigmoïde est remplacée par un softmax et la dernière couche produit $K$ sorties:

La perte est l’entropie croisée catégorielle.

Extracteur de caractéristiques et tête linéaire¶

Dans tous les cas, on peut écrire le réseau comme la composition de deux parties. Posons $\boldsymbol{\phi}(\mathbf{x}; \boldsymbol{\theta}_\phi) = \varphi(W_2\, \varphi(W_1 \mathbf{x} + \mathbf{b}_1) + \mathbf{b}_2) \in \mathbb{R}^h$, la représentation apprise par les couches cachées. Les prédictions s’écrivent alors:

C’est exactement l’équation (7) du début du chapitre. La dernière couche est un modèle linéaire (chapitre 2) ou une régression logistique (chapitre 3) appliqué aux caractéristiques apprises $\boldsymbol{\phi}(\mathbf{x})$. Les modèles des chapitres 2 et 3 sont le cas particulier $\boldsymbol{\phi}(\mathbf{x}) = \mathbf{x}$ (aucune couche cachée).

Limites du MLP¶

Le MLP traite son entrée comme un vecteur plat $\mathbf{x} \in \mathbb{R}^d$: chaque multiplication $W_\ell \mathbf{z}_{\ell-1}$ opère sur toutes les composantes de $\mathbf{z}_{\ell-1}$ sans distinction. La matrice $W_\ell$ est pleine (dense), ce qui signifie que chaque composante de la sortie dépend de toutes les composantes de l’entrée. Il n’y a aucune notion de structure spatiale ou temporelle. Pour des données tabulaires (âge, revenu, nombre de pièces), c’est approprié: il n’y a pas d’ordre naturel entre les variables.

Mais pour une image de $28 \times 28$ pixels, le MLP la transforme en un vecteur de 784 entrées. La matrice $W_1 \in \mathbb{R}^{h \times 784}$ mélange toutes les positions spatiales: elle ne sait pas que le pixel $(0, 0)$ est voisin du pixel $(0, 1)$ mais éloigné du pixel $(27, 27)$. Pour une phrase de 10 mots, le MLP a besoin d’une entrée $\mathbf{x} \in \mathbb{R}^{10d}$ de taille fixe. Que fait-on avec une phrase de 20 mots?

Ces limitations motivent des architectures qui exploitent la structure des données. Les réseaux convolutifs remplacent la matrice dense par une opération de convolution qui respecte la structure spatiale des images. Les réseaux récurrents traitent les séquences élément par élément en maintenant un état interne. Le mécanisme d’attention et les transformeurs permettent à chaque position d’une séquence de consulter directement toutes les autres, sans contrainte de longueur fixe.

Résumé¶

Ce chapitre a montré comment les réseaux de neurones s’inscrivent dans la progression des modèles vus dans les chapitres précédents. Le point de départ est toujours le cadre de maximum de vraisemblance: un modèle prédit les paramètres d’une distribution conditionnelle. La nouveauté est que la transformation des entrées (la fonction $\boldsymbol{\phi}$) est désormais apprise plutôt que fixée à l’avance. Le problème XOR a illustré pourquoi cette flexibilité est nécessaire: certaines fonctions simples sont inaccessibles aux modèles linéaires, et une couche cachée suffit à les résoudre en transformant l’espace des entrées.

La dérivation automatique calcule les gradients en décomposant un programme en opérations élémentaires et en appliquant la règle de la chaîne. Le mode arrière (VJP) produit le gradient par rapport à tous les paramètres en une seule passe, ce qui en fait la base de l’entraînement des réseaux. Les bibliothèques modernes (JAX, PyTorch) implémentent ce mécanisme automatiquement via le traçage d’opérations.

En pratique, un MLP se décompose en un extracteur de caractéristiques (les couches cachées) et une tête linéaire (la couche de sortie), ce qui généralise directement les modèles des chapitres 2 et 3. Cependant, le MLP traite ses entrées comme des vecteurs plats, sans exploiter la structure spatiale ou temporelle des données. Le chapitre suivant couvre les algorithmes d’optimisation, la stabilisation de l’entraînement et la régularisation. Les chapitres sur les réseaux récurrents, l’attention et les transformeurs présentent ensuite des architectures qui remédient aux limitations du MLP.

Exercices¶

Les exercices ★ vérifient la compréhension de base. Les exercices ★★ demandent d’appliquer les concepts à des calculs concrets. Les exercices ★★★ approfondissent le sujet et sont optionnels pour IFT3395.