Reliez l'optimisation à l'apprentissage

Quand on dit qu’un modèle « apprend », que se passe-t-il concrètement dans la machine ? Presque toujours la même chose : un algorithme d’optimisation ajuste des paramètres pour rendre une erreur la plus petite possible. La régression linéaire, la régression logistique, les SVM et les réseaux de neurones ne diffèrent que par la fonction à minimiser et par la méthode employée pour le faire.

Ce dernier chapitre relie tout ce que vous avez vu à l’apprentissage automatique : la convexité, la descente de gradient et son pas, Newton et BFGS, les multiplicateurs de Lagrange et KKT. Vous y découvrirez aussi les outils propres au domaine : le gradient stochastique, le momentum et Adam. Pour les modèles eux-mêmes, reportez-vous au cours d’Apprentissage automatique.

Entraîner un modèle, c’est minimiser

Un modèle $f_\theta$ dépend de paramètres $\theta$ (les poids d’une régression, ceux d’un réseau de neurones). On dispose de $n$ exemples $(x_i, y_i)$ et d’une fonction de perte $\ell$ qui mesure l’écart entre une prédiction et la vraie valeur. Entraîner le modèle, c’est résoudre

$$\min_\theta \; J(\theta) = \frac{1}{n} \sum_{i=1}^{n} \ell\big(f_\theta(x_i), y_i\big)$$

• $J(\theta)$ est le risque empirique, la perte moyenne sur les données d’apprentissage ;
• $\ell$ peut être le carré de l’erreur (régression), l’entropie croisée (classification), la perte charnière (SVM)…

C’est un problème sans contrainte (partie 2), souvent en très grande dimension, et dont la fonction objectif est une somme sur les exemples. Cette structure de somme va tout changer.

Le cas convexe par excellence : les moindres carrés

Pour une régression linéaire, on range les exemples dans une matrice $X$ (une ligne par exemple, une colonne par variable, plus une colonne de 1 pour la constante) et l’on minimise la somme des carrés des erreurs :

$$\min_w \; \lVert X w - y \rVert^2$$

Son gradient vaut $2X^T(Xw - y)$ et sa Hessienne $2X^T X$, semi-définie positive : le problème est convexe, c’est une quadratique du chapitre 5. Annuler le gradient donne les équations normales :

$$X^T X \, w = X^T y$$

Si les colonnes de $X$ sont linéairement indépendantes, $X^T X$ est définie positive et la solution est unique.

Exemple

Ajustons une droite $y = a x + b$ aux trois points $(0, 1)$, $(1, 3)$ et $(2, 4)$. Avec les sommes $\sum x_i = 3$, $\sum x_i^2 = 5$, $\sum y_i = 8$ et $\sum x_i y_i = 11$, les équations normales s’écrivent

$$5a + 3b = 11, \qquad 3a + 3b = 8$$

En soustrayant, $2a = 3$, donc $a = 1{,}5$ et $b = 7/6 \approx 1{,}167$. Les résidus valent $\tfrac16$, $-\tfrac13$ et $\tfrac16$ : leur somme est nulle, comme toujours quand le modèle contient une constante.

La droite des moindres carrés rend minimale la somme des carrés des segments verticaux (les résidus) : 1/6, 1/3 et 1/6.

Attention

En pratique, on ne forme pas $X^T X$ pour la résoudre : son conditionnement est le carré de celui de $X$, ce qui amplifie les erreurs d’arrondi. On utilise une factorisation QR ou np.linalg.lstsq. Et quand $n$ ou le nombre de variables devient énorme, on revient… à la descente de gradient.

Le gradient stochastique

Le gradient de $J$ est une moyenne de $n$ gradients individuels : $\nabla J(\theta) = \frac1n \sum_i \nabla \ell_i(\theta)$. Avec un million d’exemples, chaque pas de descente de gradient demande de parcourir tout le jeu de données. C’est beaucoup de calcul pour un seul petit pas.

L’idée du gradient stochastique (stochastic gradient descent, SGD) : estimer ce gradient sur un petit échantillon tiré au hasard, un mini-lot (mini-batch) $B$ de quelques dizaines d’exemples :

$$\theta_{k+1} = \theta_k - \alpha_k \, \frac{1}{\lvert B_k \rvert} \sum_{i \in B_k} \nabla \ell_i(\theta_k)$$

Cette estimation est sans biais (en moyenne, elle vaut le vrai gradient) mais bruitée. On fait donc beaucoup plus de pas, chacun beaucoup moins cher, et l’on parcourt le jeu de données par époques (une époque = tous les exemples vus une fois).

	Gradient complet	SGD (un exemple)	Mini-lots
Coût d’un pas	$n$ gradients	1 gradient	$\lvert B \rvert$ gradients (32 à 512)
Bruit	aucun	très fort	modéré
Convergence avec pas fixe	vers le minimum	fluctue autour du minimum	fluctue, moins fort
Usage	petits jeux de données	historique	standard en apprentissage profond

Le bruit a deux conséquences. D’une part, avec un pas fixe, les itérées ne convergent pas exactement : elles tournent autour du minimum dans une zone d’autant plus large que le pas est grand. Pour converger, il faut un pas décroissant (par exemple $\alpha_k \propto 1/k$). D’autre part, ce bruit aide à s’échapper des points selles, très nombreux dans les réseaux de neurones, où un gradient exact resterait coincé.

Sur la régression du code ci-dessous : le SGD suit en gros la même route que le gradient complet, mais en tremblant. Avec un pas fixe, il finit par fluctuer autour du minimum ; un pas décroissant calme ces fluctuations.

Accélérer : momentum et Adam

Le momentum

Au chapitre 6, la descente de gradient zigzaguait dans les vallées étroites. Le momentum (méthode de la boule lourde) lisse la trajectoire en accumulant une « vitesse » :

$$v_{k+1} = \beta \, v_k + \nabla J(\theta_k), \qquad \theta_{k+1} = \theta_k - \alpha \, v_{k+1}$$

avec $\beta$ proche de 1 (typiquement 0,9). Dans les directions où le gradient change de signe à chaque pas (le zigzag), les contributions successives se compensent ; dans la direction de la vallée, elles s’additionnent. La boule prend de l’élan dans la bonne direction.

Dans une vallée étroite, le gradient seul rebondit d'un bord à l'autre. Avec le momentum, les allers-retours se compensent et l'élan s'accumule le long de la vallée : en 40 itérations, il arrive au fond.

Adam

Adam (adaptive moment estimation) combine le momentum et un pas adapté à chaque coordonnée. Avec $g_k$ le gradient (stochastique) à l’itération $k$ :

$$m_k = \beta_1 m_{k-1} + (1 - \beta_1) g_k, \qquad s_k = \beta_2 s_{k-1} + (1 - \beta_2) g_k^2$$ $$\theta_{k+1} = \theta_k - \alpha \, \frac{\hat m_k}{\sqrt{\hat s_k} + \varepsilon}, \qquad \hat m_k = \frac{m_k}{1 - \beta_1^k}, \quad \hat s_k = \frac{s_k}{1 - \beta_2^k}$$

• $m_k$ est une moyenne mobile des gradients (le momentum) ;
• $s_k$ est une moyenne mobile de leurs carrés, coordonnée par coordonnée ;
• $\hat m_k$ et $\hat s_k$ corrigent le biais dû à l’initialisation à zéro ;
• $\varepsilon$ (de l’ordre de $10^{-8}$) évite la division par zéro.

Diviser par $\sqrt{\hat s_k}$ donne de grands pas aux paramètres dont le gradient est faible et de petits pas à ceux dont le gradient est fort. C’est une forme bon marché de préconditionnement : Adam compense en partie le mauvais conditionnement, un peu comme Newton, mais sans Hessienne. Les valeurs par défaut ($\alpha = 10^{-3}$, $\beta_1 = 0{,}9$, $\beta_2 = 0{,}999$) fonctionnent souvent sans réglage, ce qui en a fait l’optimiseur par défaut de l’apprentissage profond.

À retenir

Le taux d’apprentissage joue le rôle du pas du chapitre 6 : trop grand, l’entraînement oscille ou diverge ; trop petit, il stagne. Le SGD remplace le gradient par une estimation sur un mini-lot ; le momentum lisse la trajectoire ; Adam adapte en plus le pas à chaque paramètre.

Voici un SGD par mini-lots avec momentum sur une régression linéaire, comparé à la solution exacte des moindres carrés :

import numpy as np

rng = np.random.default_rng(0)
n = 1000
X = np.column_stack([np.ones(n), rng.uniform(0, 10, n)])   # constante + une variable
y = 2.0 + 0.5 * X[:, 1] + rng.normal(0, 1, n)               # vraie droite : 2 + 0,5 x

w_exact, *_ = np.linalg.lstsq(X, y, rcond=None)              # moindres carrés exacts

w, v = np.zeros(2), np.zeros(2)
pas, beta, taille_lot = 0.001, 0.9, 32
for epoque in range(50):
    ordre = rng.permutation(n)                               # on mélange à chaque époque
    for debut in range(0, n, taille_lot):
        lot = ordre[debut:debut + taille_lot]
        grad = 2 * X[lot].T @ (X[lot] @ w - y[lot]) / len(lot)   # gradient sur le mini-lot
        v = beta * v + grad                                      # momentum
        w = w - pas * v

print(w_exact.round(3), w.round(3))   # [1.999 0.498] [1.985 0.474]
# les deux sont proches de (2, 0.5) ; avec un pas fixe, w fluctue de quelques centièmes d'une époque à l'autre

Avec PyTorch, la même boucle se réduit à quelques lignes, l’optimiseur se chargeant des mises à jour :

import torch

modele = torch.nn.Linear(1, 1)
optimiseur = torch.optim.Adam(modele.parameters(), lr=1e-2)
perte = torch.nn.MSELoss()

x = torch.rand(1000, 1) * 10
y = 2 + 0.5 * x + torch.randn(1000, 1)
for epoque in range(200):
    optimiseur.zero_grad()          # remet les gradients à zéro
    erreur = perte(modele(x), y)    # risque empirique
    erreur.backward()               # calcule le gradient (rétropropagation)
    optimiseur.step()               # un pas d'Adam

La régularisation, une contrainte déguisée

Pour éviter le surapprentissage, on régularise : on empêche les poids de devenir trop grands. Deux formulations semblent différentes.

• Contrainte : $\min_w J(w)$ sous $\lVert w \rVert^2 \leq t$ (un budget sur la taille des poids).
• Pénalité : $\min_w J(w) + \lambda \lVert w \rVert^2$ (la régression Ridge quand $J$ est l’erreur quadratique).

Les conditions KKT du chapitre 9 montrent qu’elles sont équivalentes. Pour la version contrainte, la stationnarité s’écrit $\nabla J(w^*) + 2\mu w^* = 0$, qui est exactement la condition d’optimalité de la version pénalisée avec $\lambda = \mu$. Le coefficient de régularisation est un multiplicateur de Lagrange. À chaque budget $t$ correspond un $\lambda$, et inversement : un grand $\lambda$ correspond à un petit budget.

Pour les moindres carrés, la solution Ridge vérifie

$$(X^T X + \lambda I) \, w = X^T y$$

Ajouter $\lambda I$ augmente toutes les valeurs propres de $\lambda$ : la matrice devient définie positive même si les variables sont corrélées, et son conditionnement s’améliore. La régularisation stabilise donc à la fois le modèle et l’optimisation.

Avec la norme $\ell_1$ ($\lVert w \rVert_1 = \sum_j \lvert w_j \rvert \leq t$), l’ensemble réalisable est un losange (en dimension 2) dont les coins sont sur les axes. Les lignes de niveau de $J$ le touchent souvent en un coin, où certains poids sont exactement nuls : c’est pourquoi le Lasso sélectionne des variables.

Mêmes lignes de niveau, deux budgets. Le cercle n'a pas de coin : le contact tombe en général hors des axes. Le losange a des coins sur les axes : le contact y tombe souvent, et un poids s'annule exactement.

Piège

La régularisation pénalise la taille des poids, donc elle dépend de l’échelle des variables. Sans standardisation préalable, une variable exprimée en euros (poids minuscule) serait à peine pénalisée, et une variable en milliers d’euros beaucoup plus : le résultat dépendrait d’un simple choix d’unité.

Convexité, SVM et conditions KKT

La convexité (chapitre 5) sépare les modèles en deux familles.

Modèle	Problème	Conséquence	Optimiseur typique
Régression linéaire, Ridge	quadratique convexe	minimum global, formule fermée possible	équations normales, QR, gradient
Régression logistique	convexe	minimum global, indépendant de l’initialisation	L-BFGS, Newton
SVM	quadratique convexe sous contraintes	minimum global, dual via KKT	solveurs dédiés (SMO)
Réseau de neurones	non convexe	minimums locaux, points selles	SGD, momentum, Adam

Pour un réseau de neurones, on ne vise pas le minimum global : on cherche un bon minimum local, qui généralise bien. Le résultat dépend alors de l’initialisation, du taux d’apprentissage et même du bruit du SGD.

Les SVM sont un bel exemple de KKT à l’œuvre. Pour des données séparables, la SVM à marge maximale résout

$$\min_{w, b} \; \tfrac12 \lVert w \rVert^2 \quad \text{s.c.} \quad y_i (w^T x_i + b) \geq 1 \quad \text{pour tout } i$$

où $y_i \in \lbrace -1, +1 \rbrace$ est la classe de l’exemple $x_i$. Avec un multiplicateur $\alpha_i \geq 0$ par contrainte, les conditions KKT donnent :

• stationnarité : $w = \sum_i \alpha_i y_i x_i$ et $\sum_i \alpha_i y_i = 0$ ;
• complémentarité : $\alpha_i \big(y_i (w^T x_i + b) - 1\big) = 0$.

La complémentarité a une conséquence remarquable : seuls les exemples situés exactement sur la marge peuvent avoir $\alpha_i \gt 0$. Ce sont les vecteurs de support, et $w$ ne dépend que d’eux ; tous les autres exemples pourraient disparaître sans changer le modèle. Enfin, le problème dual (chapitre 9) ne fait intervenir les données que par leurs produits scalaires $x_i^T x_j$, qu’on peut remplacer par un noyau $K(x_i, x_j)$ : c’est l’astuce du noyau.

Seuls les exemples posés exactement sur la marge (entourés) ont un multiplicateur αᵢ > 0 : ce sont les vecteurs de support. On pourrait supprimer tous les autres sans déplacer la frontière.

En entretien

« Entraîner un modèle, c’est minimiser un risque empirique. Comme c’est une somme sur les exemples, on estime le gradient sur des mini-lots : c’est le SGD, que momentum et Adam accélèrent. La régularisation Ridge est une contrainte sur la norme des poids, dont le coefficient est le multiplicateur de Lagrange. »

En résumé

• Entraîner un modèle revient à minimiser un risque empirique, une moyenne de pertes sur les exemples.
• Les moindres carrés sont convexes ; leur solution vérifie les équations normales $X^T X w = X^T y$, qu’on résout en pratique par QR ou lstsq.
• Le SGD estime le gradient sur un mini-lot : des pas bien moins chers mais bruités, qui demandent un pas décroissant pour converger exactement.
• Le momentum lisse le zigzag ; Adam ajoute un pas adapté à chaque paramètre, une forme de préconditionnement.
• La régularisation Ridge équivaut à une contrainte sur la norme des poids ; son coefficient est un multiplicateur de Lagrange, et elle améliore le conditionnement.
• Dans une SVM, la complémentarité KKT fait apparaître les vecteurs de support, et le dual permet l’astuce du noyau.

Exercices

Exercice · Résoudre les équations normales

facile

Ajustez une droite $y = ax + b$ aux points $(0, 0)$, $(1, 1)$ et $(2, 3)$ par les moindres carrés.

Voir la solution

Les sommes valent $\sum x_i = 3$, $\sum x_i^2 = 5$, $\sum y_i = 4$ et $\sum x_i y_i = 0 + 1 + 6 = 7$, pour $n = 3$. Les équations normales s’écrivent $5a + 3b = 7$ et $3a + 3b = 4$. En soustrayant, $2a = 3$, donc $a = 1{,}5$, puis $b = (4 - 4{,}5)/3 = -\tfrac16$.

La droite est $y = 1{,}5x - \tfrac16$. Résidus : $0 - (-\tfrac16) = \tfrac16$, $1 - \tfrac43 = -\tfrac13$, $3 - \tfrac{17}{6} = \tfrac16$, de somme nulle.

En code

Python

import numpy as np

x = np.array([0, 1, 2]); y = np.array([0, 1, 3])
X = np.column_stack([x, np.ones(3)])            # colonnes : a (pente), b (constante)
print(X.T @ X, X.T @ y)                         # [[5 3] [3 3]] [7 4] : équations normales
w, *_ = np.linalg.lstsq(X, y, rcond=None)
print(w)                                        # [ 1.5 -0.1667] : y = 1,5 x - 1/6
residus = y - X @ w
print(residus, residus.sum())                   # [ 0.1667 -0.3333  0.1667] 0 (à 1e-15 près)

R

x <- c(0, 1, 2); y <- c(0, 1, 3)
X <- cbind(a = x, b = 1)                        # colonnes : pente, constante
crossprod(X); crossprod(X, y)                   # [5 3; 3 3] et (7, 4) : équations normales
solve(crossprod(X), crossprod(X, y))            # a = 1.5, b = -0.1667
modele <- lm(y ~ x)                             # même résultat avec lm
coef(modele); residuals(modele)                 # -0.1667 1.5 ; 0.1667 -0.3333 0.1667

Exercice · Ridge et conditionnement

moyen

Dans une régression, la matrice $X^T X$ a pour valeurs propres 1 et 100. On la résout par descente de gradient avec le meilleur pas fixe.

1. Quel est le facteur de réduction de l’erreur par itération ?
2. Même question pour Ridge avec $\lambda = 10$.
3. Combien d’itérations faut-il, environ, pour diviser l’erreur par 1 000 dans chaque cas ?

Voir la solution

1. Le conditionnement vaut $\kappa = 100$, d’où un facteur $\frac{\kappa - 1}{\kappa + 1} = \frac{99}{101} \approx 0{,}980$.
2. Avec Ridge, la matrice devient $X^T X + 10 I$, de valeurs propres 11 et 110 : $\kappa = 10$, et le facteur vaut $\frac{9}{11} \approx 0{,}818$.
3. Il faut $k$ tel que $0{,}980^k \leq 10^{-3}$, soit $k \geq \ln(10^{-3}) / \ln(0{,}980) \approx 345$ itérations sans régularisation, contre $\ln(10^{-3}) / \ln(0{,}818) \approx 35$ avec Ridge. La régularisation rend l’optimisation environ dix fois plus rapide (au prix d’un léger biais sur les coefficients).

En code

Python

import numpy as np

for nom, vp in [("sans régularisation", [1, 100]), ("Ridge λ = 10", [11, 110])]:
    lmin, lmax = min(vp), max(vp)
    facteur = (lmax / lmin - 1) / (lmax / lmin + 1)
    # Vérification : gradient au meilleur pas fixe sur la quadratique ½ w'Aw
    A, pas = np.diag(vp).astype(float), 2 / (lmin + lmax)
    w, k = np.ones(2), 0
    while np.linalg.norm(w) > 1e-3 * np.sqrt(2):   # erreur divisée par 1000
        w, k = w - pas * A @ w, k + 1
    print(nom, round(facteur, 3), round(np.log(1e-3) / np.log(facteur), 1), k)
# sans régularisation 0.98 345.4 346
# Ridge λ = 10 0.818 34.4 35

R

for (vp in list(c(1, 100), c(11, 110))) {
  kappa <- max(vp) / min(vp)
  facteur <- (kappa - 1) / (kappa + 1)
  # Vérification : gradient au meilleur pas fixe sur la quadratique ½ w'Aw
  A <- diag(vp); pas <- 2 / sum(vp)
  w <- c(1, 1); k <- 0
  while (sqrt(sum(w^2)) > 1e-3 * sqrt(2)) { w <- w - pas * A %*% w; k <- k + 1 }
  cat("kappa", kappa, ": facteur", round(facteur, 3),
      ", théorie", round(log(1e-3) / log(facteur), 1), ", simulation", k, "\n")
}
# kappa 100 : facteur 0.98 , théorie 345.4 , simulation 346
# kappa 10 : facteur 0.818 , théorie 34.4 , simulation 35

Exercice · Trouver les vecteurs de support

difficile

On dispose de deux exemples, $x_1 = (1, 1)$ de classe $+1$ et $x_2 = (-1, -1)$ de classe $-1$.

1. Résolvez la SVM à marge maximale (on admettra que $b = 0$ par symétrie). Quelle est la largeur de la marge $2 / \lVert w \rVert$ ?
2. Calculez les multiplicateurs $\alpha_1$ et $\alpha_2$ avec les conditions KKT.
3. On ajoute un troisième exemple $x_3 = (2, 3)$ de classe $+1$. Que deviennent la solution et $\alpha_3$ ?

Voir la solution

1. Par symétrie, $w = (a, a)$ avec $b = 0$. Les contraintes s’écrivent $w^T x_1 = 2a \geq 1$ et $-w^T x_2 = 2a \geq 1$. On minimise $\tfrac12 \lVert w \rVert^2 = a^2$, donc on prend $a$ le plus petit possible : $a = \tfrac12$, et $w = (0{,}5 ; 0{,}5)$. La marge vaut $2 / \lVert w \rVert = 2 / 0{,}707 \approx 2{,}83$, exactement la distance entre les deux points ($2\sqrt{2}$) : la frontière passe au milieu.
2. La stationnarité donne $\alpha_1 = \alpha_2 = \alpha$ (car $\sum \alpha_i y_i = 0$) et $w = \alpha x_1 - \alpha x_2 = 2\alpha (1, 1)$. Comme $w = (0{,}5 ; 0{,}5)$, $\alpha = \tfrac14$. Les deux multiplicateurs sont strictement positifs : les deux points sont des vecteurs de support, et les deux contraintes sont bien actives.
3. Pour $x_3$, $w^T x_3 = 0{,}5 \times 2 + 0{,}5 \times 3 = 2{,}5 \gt 1$ : la contrainte est inactive. Par complémentarité, $\alpha_3 = 0$. Les conditions KKT restent vérifiées avec les mêmes $w$, $b$, $\alpha_1$ et $\alpha_2$, et le problème est convexe : la solution ne change pas. Le point $x_3$, loin de la frontière, n’est pas un vecteur de support.

En code

Python

import numpy as np
from sklearn.svm import SVC

X = np.array([[1, 1], [-1, -1]]); y = np.array([1, -1])
svm = SVC(kernel="linear", C=1e6).fit(X, y)     # C très grand : marge dure
w = svm.coef_[0]
print(w, svm.intercept_, 2 / np.linalg.norm(w))  # [0.5 0.5] [0.] 2.83
print(svm.dual_coef_)                            # [[-0.25  0.25]] : α_i y_i, donc α = 1/4

# On ajoute x3 = (2, 3) de classe +1
X3 = np.vstack([X, [2, 3]]); y3 = np.append(y, 1)
svm3 = SVC(kernel="linear", C=1e6).fit(X3, y3)
print(svm3.coef_[0], svm3.support_)              # [0.5 0.5] [1 0] : x3 n'est pas support
print(w @ [2, 3])                                # 2.5 > 1 : contrainte inactive, α3 = 0

R

X <- rbind(c(1, 1), c(-1, -1), c(2, 3)); y <- c(1, -1, 1)
svm <- function(n) {                  # marge dure sur les n premiers points
  ui <- y[1:n] * cbind(X[1:n, , drop = FALSE], 1)   # y_i (w'x_i + b) - 1 >= 0
  res <- constrOptim(c(2, 2, 0), function(t) sum(t[1:2]^2) / 2,
                     grad = function(t) c(t[1:2], 0), ui = ui, ci = rep(1, n),
                     outer.eps = 1e-10)
  round(res$par, 3)                   # (w1, w2, b)
}
t2 <- svm(2); t2                      # 0.5 0.5 0 : w = (0.5, 0.5), b = 0
2 / sqrt(sum(t2[1:2]^2))              # 2.83 : largeur de la marge
# Stationnarité w = α (x1 - x2) = 2α (1, 1), donc α1 = α2 = 1/4
t2[1] / 2
svm(3)                                # 0.5 0.5 0 : la solution ne change pas
sum(t2[1:2] * X[3, ])                 # 2.5 > 1 : x3 hors de la marge, α3 = 0