Accélérez avec Newton et quasi-Newton

La descente de gradient n’utilise qu’une information : la pente. Elle ignore si le terrain est une cuvette étroite ou une large plaine, et c’est pour cela qu’elle zigzague dès que les lignes de niveau s’allongent. Notre randonneur ferait mieux s’il sentait aussi la courbure du sol : dans une cuvette bien arrondie, il pourrait deviner où se trouve le fond et y sauter directement.

C’est l’idée de la méthode de Newton. Elle converge en une poignée d’itérations là où le gradient en demande des centaines, mais elle coûte plus cher et peut s’égarer. Dans ce chapitre, vous allez la voir à l’œuvre sur la fonction de notre TP, comprendre ses limites, puis découvrir les méthodes de quasi-Newton, le compromis utilisé en pratique.

L’idée : sauter au fond d’un modèle quadratique

Près du point courant $x_k$, Taylor à l’ordre 2 approche $f$ par une quadratique :

$$f(x_k + h) \approx f(x_k) + \nabla f(x_k)^T h + \tfrac12 \, h^T \nabla^2 f(x_k) \, h$$

Si la Hessienne $\nabla^2 f(x_k)$ est définie positive, ce modèle a un unique minimum, là où son gradient $\nabla f(x_k) + \nabla^2 f(x_k) h$ s’annule. La méthode de Newton saute directement à ce minimum :

$$x_{k+1} = x_k - \big[\nabla^2 f(x_k)\big]^{-1} \nabla f(x_k)$$

• $\nabla f(x_k)$ est le gradient au point courant ;
• $\nabla^2 f(x_k)$ est la Hessienne, et $[\cdot]^{-1}$ son inverse ;
• le déplacement $h_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$ s’appelle la direction de Newton.

En dimension 1, la formule devient $x_{k+1} = x_k - f'(x_k) / f''(x_k)$. Vous reconnaissez peut-être la méthode de Newton pour résoudre une équation : c’est exactement elle, appliquée à l’équation $f'(x) = 0$. On trace la tangente au graphe de $f'$ et l’on prend son intersection avec l’axe horizontal.

Une itération de Newton : on remplace f par sa parabole de Taylor en x₀ = 1 et l'on saute au fond de cette parabole, x₁ = 1 − f′(1)/f″(1) = 1,25.

Comparez avec le gradient : Newton, c’est un gradient dont le pas est remplacé par l’inverse de la courbure. Là où la fonction est très courbée, il fait de petits pas ; là où elle est plate, de grands pas. Et comme il corrige chaque direction selon sa propre courbure, le conditionnement ne le gêne plus.

Newton en dimension 1 : notre TP

En TP, nous avons appliqué Newton à $f(x) = x^4 - 7x + 8$, avec $f'(x) = 4x^3 - 7$ et $f''(x) = 12x^2$, depuis trois points de départ : 0,8, 1 et 1,4. Avec les mêmes critères d’arrêt que pour le gradient (tolérance 0,025 sur $\lvert f'(x) \rvert$, 20 itérations au plus), la sortie affichée pour le dernier départ, $x_0 = 1{,}4$, donne un minimum en $x \approx 1{,}205611$, avec $f(x) \approx 1{,}673379$. Si l’on rejoue le code du TP avec ces critères, Newton s’arrête après deux itérations seulement depuis 1,4 (et trois depuis 0,8 ou 1), là où le gradient à pas fixe 0,1 en enchaînait 18 pour arriver au même endroit.

Pour voir la vitesse de Newton, suivons l’erreur $\lvert x_k - x^* \rvert$ depuis $x_0 = 1$, en continuant au-delà de la tolérance du TP ($x^* = (7/4)^{1/3} \approx 1{,}205071$) :

$k$	$x_k$	Erreur
0	1	$2{,}1 \times 10^{-1}$
1	1,25	$4{,}5 \times 10^{-2}$
2	1,206667	$1{,}6 \times 10^{-3}$
3	1,205073	$2{,}1 \times 10^{-6}$
4	1,205071	$3{,}7 \times 10^{-12}$

Le nombre de décimales exactes double à chaque itération. C’est la convergence quadratique : près de la solution, l’erreur vérifie $\lvert x_{k+1} - x^* \rvert \leq C \lvert x_k - x^* \rvert^2$ pour une constante $C$. En comparaison, le gradient a une convergence linéaire : l’erreur est multipliée par un facteur constant à chaque pas, et chaque décimale coûte le même nombre d’itérations.

En échelle logarithmique, la convergence linéaire du gradient est une droite ; celle de Newton est une chute de plus en plus raide, car le nombre de décimales exactes double à chaque pas.

À retenir

La méthode de Newton minimise à chaque itération le modèle quadratique de Taylor : $x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$. Près d’un minimum à Hessienne définie positive, elle converge quadratiquement : le nombre de chiffres exacts double à chaque pas.

Quand Newton part trop loin

Ce bel ordre ne vaut que près de la solution. Partons de $x_0 = 0{,}1$ : la courbure y est minuscule, $f''(0{,}1) = 0{,}12$, alors que la pente est forte, $f'(0{,}1) \approx -7$. Le modèle quadratique est une parabole très évasée dont le fond est loin, et Newton y saute : $x_1 = 0{,}1 + 6{,}996 / 0{,}12 \approx 58{,}4$. Loin de l’origine, $f \approx x^4$ et chaque pas de Newton ne réduit $x$ que d’un tiers : 38,9, puis 26,0, puis 17,3… Il faut une douzaine d’itérations pour revenir près de 1,2. Et depuis $x_0 = 0$, c’est pire : $f''(0) = 0$, la formule divise par zéro et la méthode ne peut même pas démarrer.

Parti de 0,1, où la courbure est minuscule, Newton saute à 58,4 puis ne revient que d'un tiers à chaque pas : il faut une douzaine d'itérations avant que la convergence rapide ne démarre.

Newton en dimension n

En dimension $n$, on ne calcule jamais l’inverse de la Hessienne. On résout le système linéaire

$$\nabla^2 f(x_k) \, h_k = -\nabla f(x_k)$$

puis on pose $x_{k+1} = x_k + h_k$. Résoudre un système coûte moins cher et est numériquement plus stable qu’inverser une matrice.

Exemple

Appliquons Newton à la quadratique $f_1(x, y) = 2x^2 + 3x + y^2 - 2$ du TP, depuis $(10, 10)$. Le gradient vaut $(43, 20)$ et la Hessienne $\mathrm{diag}(4, 2)$. Le système donne $h = (-43/4, -20/2) = (-10{,}75 ; -10)$, et le nouveau point est $(-0{,}75 ; 0)$ : le minimum exact, en une seule itération.

Ce n’est pas un hasard : pour une quadratique, le modèle de Taylor est la fonction elle-même, et Newton tombe directement sur son minimum. Le gradient à pas optimal, lui, avait zigzagué pendant plusieurs itérations, et le pas fixe 0,1 en demandait plusieurs dizaines pour atteindre une grande précision.

Sur une quadratique, Newton corrige chaque direction selon sa courbure et tombe sur le minimum en une itération, là où le gradient zigzague.

Le revers de la médaille, c’est le coût. À chaque itération, il faut :

• calculer la Hessienne, soit $n^2$ dérivées secondes (environ $n^2 / 2$ grâce à la symétrie) ;
• résoudre un système $n \times n$, en $O(n^3)$ opérations.

Pour $n = 100$, c’est instantané. Pour un réseau de neurones à $n = 10^7$ paramètres, la Hessienne aurait $10^{14}$ coefficients : impossible même à stocker.

Les limites de Newton, et comment les contourner

Newton cherche un point où le gradient s’annule. Il ne sait pas faire la différence entre un minimum, un maximum et un col. Si la Hessienne n’est pas définie positive, la direction de Newton n’est pas forcément une direction de descente, et la méthode peut monter.

Piège

Reprenons l’altitude du plateau, $g(x, y) = x^3 - 3x + y^2$, qui a un col en $(-1, 0)$. Depuis $(-0{,}9 ; 0)$, le gradient vaut $(-0{,}57 ; 0)$ et la Hessienne $\mathrm{diag}(-5{,}4 ; 2)$. Le pas de Newton donne $x_1 = (0{,}81 + 1) / (-1{,}8) \approx -1{,}0056$ : l’algorithme se précipite vers le col, et l’altitude augmente, de 1,971 à 1,9999. Newton a trouvé un point critique… qui n’est pas un minimum.

Newton vise le point où le gradient s'annule, que ce soit un minimum ou un col : ici, il remonte vers le col, alors que le gradient indique la descente vers l'est.

Deux remèdes classiques rendent la méthode robuste :

• Newton amorti (damped Newton) : on garde la direction de Newton mais l’on choisit le pas par une recherche linéaire d’Armijo (chapitre 6). Loin de la solution, on fait des pas raisonnables ; près d’elle, le pas 1 est accepté et l’on retrouve la convergence quadratique.
• Correction de la Hessienne : si elle n’est pas définie positive, on la remplace par $\nabla^2 f(x_k) + \lambda I$ avec $\lambda \gt 0$ assez grand. Quand $\lambda$ est grand, on retombe sur un petit pas de gradient ; quand il est petit, sur Newton. C’est l’idée de la méthode de Levenberg-Marquardt, très utilisée pour les moindres carrés non linéaires.

Les méthodes de quasi-Newton

Peut-on avoir la vitesse de Newton sans payer la Hessienne ? Presque. Les méthodes de quasi-Newton construisent une approximation de la Hessienne (ou de son inverse) à partir des gradients déjà calculés, sans aucune dérivée seconde.

L’idée vient de la dimension 1 : la dérivée seconde est la variation de la dérivée, que l’on peut estimer entre deux itérées. En dimension $n$, on note $s_k = x_{k+1} - x_k$ le déplacement et $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$ la variation du gradient. On impose à la nouvelle approximation $B_{k+1}$ de la Hessienne la condition de la sécante :

$$B_{k+1} \, s_k = y_k$$

Autrement dit, $B_{k+1}$ doit expliquer la variation de gradient qu’on vient d’observer. Parmi toutes les matrices qui la vérifient, la méthode BFGS (Broyden, Fletcher, Goldfarb, Shanno) choisit la plus proche de la précédente, et met directement à jour l’approximation $H_k$ de l’inverse de la Hessienne :

$$H_{k+1} = \big(I - \rho_k s_k y_k^T\big) H_k \big(I - \rho_k y_k s_k^T\big) + \rho_k s_k s_k^T, \qquad \rho_k = \frac{1}{y_k^T s_k}$$

La direction est alors $-H_k \nabla f(x_k)$ : un simple produit matrice-vecteur, en $O(n^2)$, sans système à résoudre. Si $H_k$ est définie positive et $y_k^T s_k \gt 0$ (ce que garantit une bonne recherche linéaire), $H_{k+1}$ l’est aussi : les directions restent des directions de descente. La convergence est superlinéaire : plus rapide que linéaire, un peu moins que quadratique.

Pour les très grandes dimensions, L-BFGS (limited memory) ne stocke même pas la matrice $H_k$ : il garde seulement les $m$ derniers couples $(s_k, y_k)$, avec $m$ de l’ordre de 5 à 20, et reconstruit le produit $H_k \nabla f$ à la volée. C’est l’algorithme par défaut de nombreux solveurs, notamment pour la régression logistique dans scikit-learn.

Exemple

Comparons les trois méthodes sur la fonction de Rosenbrock, $f(x, y) = (1 - x)^2 + 100 (y - x^2)^2$, un test classique : son minimum $(1, 1)$ est au fond d’une vallée étroite et courbée, très mal conditionnée. Départ $(-1{,}2 ; 1)$, arrêt quand $\lVert \nabla f \rVert \lt 10^{-6}$, avec une implémentation simple de chaque méthode (recherche d’Armijo pour le gradient et BFGS) :

Méthode	Itérations
Gradient avec Armijo	13 756
BFGS avec Armijo	34
Newton (pas 1)	6

Le gradient s’épuise à zigzaguer dans la vallée. BFGS fait presque aussi bien que Newton, sans jamais calculer de dérivée seconde.

Sur la vallée courbe de Rosenbrock, le gradient s'épuise en 13 756 itérations, BFGS suit la vallée en 34, Newton conclut en 6 malgré un grand saut hors du cadre.

En pratique, on n’implémente pas BFGS soi-même, on appelle SciPy :

import numpy as np
from scipy.optimize import minimize

def rosenbrock(v):
    x, y = v
    return (1 - x) ** 2 + 100 * (y - x ** 2) ** 2

def grad_rosenbrock(v):
    x, y = v
    return np.array([-2 * (1 - x) - 400 * x * (y - x ** 2), 200 * (y - x ** 2)])

res = minimize(rosenbrock, x0=[-1.2, 1.0], jac=grad_rosenbrock, method="BFGS")
print(res.x, res.nit)   # un point très proche de (1, 1), en quelques dizaines d'itérations

Pour choisir, retenez ce tableau :

	Gradient	Newton	Quasi-Newton (BFGS)
Information utilisée	gradient	gradient et Hessienne	gradients successifs
Coût d’une itération	$O(n)$	$O(n^3)$	$O(n^2)$, ou $O(mn)$ pour L-BFGS
Vitesse près de la solution	linéaire, sensible au conditionnement	quadratique	superlinéaire
Robustesse loin de la solution	bonne avec un petit pas	faible sans amortissement	bonne avec recherche linéaire
Usage typique	très grandes dimensions, apprentissage profond	petites dimensions, haute précision	dimensions moyennes à grandes

En entretien

« Newton utilise la Hessienne pour corriger la direction selon la courbure : il converge en très peu d’itérations, mais chaque itération coûte un système linéaire et il peut aller vers un col. BFGS approche la Hessienne à partir des gradients successifs : c’est le bon compromis en pratique, et L-BFGS le rend utilisable en grande dimension. »

En résumé

• Newton minimise le modèle quadratique de Taylor : $x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$ ; en dimension 1, $x_{k+1} = x_k - f'(x_k) / f''(x_k)$.
• Près d’un minimum, sa convergence est quadratique ; sur une quadratique, il trouve le minimum en une itération.
• On résout le système $\nabla^2 f(x_k) h_k = -\nabla f(x_k)$ sans inverser la Hessienne ; le coût est en $O(n^3)$ par itération.
• Loin de la solution, Newton peut faire des sauts énormes, et il converge vers n’importe quel point critique, col compris ; on le stabilise par recherche linéaire ou en ajoutant $\lambda I$ à la Hessienne.
• Les méthodes de quasi-Newton (BFGS, L-BFGS) approchent la Hessienne grâce à la condition de la sécante $B_{k+1} s_k = y_k$ : convergence superlinéaire sans dérivée seconde.

Exercices

Exercice · Trois itérations de Newton

facile

Appliquez la méthode de Newton à $f(x) = e^x - 2x$ depuis $x_0 = 0$ pendant trois itérations. Quel est le minimum exact ?

Voir la solution

On a $f'(x) = e^x - 2$ et $f''(x) = e^x \gt 0$ : $f$ est strictement convexe, et son minimum vérifie $e^x = 2$, soit $x^* = \ln 2 \approx 0{,}693147$.

La formule donne $x_{k+1} = x_k - (e^{x_k} - 2)/e^{x_k} = x_k - 1 + 2e^{-x_k}$.

• $x_1 = 0 - 1 + 2 = 1$ ;
• $x_2 = 1 - 1 + 2e^{-1} \approx 0{,}735759$ ;
• $x_3 \approx 0{,}694042$.

Les erreurs valent environ 0,31, puis 0,043, puis 0,0009 : le nombre de décimales exactes double à peu près à chaque pas. L’itération suivante donne 0,693148, à $10^{-6}$ près.

En code

Python

import numpy as np

x = 0.0
for k in range(1, 5):
    x = x - (np.exp(x) - 2) / np.exp(x)      # x - f'(x) / f''(x)
    print(k, round(x, 6), f"erreur {abs(x - np.log(2)):.1e}")
# 1 1.0      erreur 3.1e-01
# 2 0.735759 erreur 4.3e-02
# 3 0.694042 erreur 9.0e-04
# 4 0.693148 erreur 4.0e-07

R

x <- 0
for (k in 1:4) {
  x <- x - (exp(x) - 2) / exp(x)             # x - f'(x) / f''(x)
  cat(k, round(x, 6), "erreur", signif(abs(x - log(2)), 2), "\n")
}
# 1 1 erreur 0.31
# 2 0.735759 erreur 0.043
# 3 0.694042 erreur 9e-04
# 4 0.693148 erreur 4e-07

Exercice · Newton sur une quadratique

moyen

Appliquez une itération de Newton à $f(x, y) = x^2 + xy + y^2 - 3x$ depuis l’origine. Comparez au minimum trouvé au chapitre 5, et expliquez.

Voir la solution

Le gradient vaut $(2x + y - 3, \; x + 2y)$, soit $(-3, 0)$ à l’origine. La Hessienne est constante : $H = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$. On résout $H h = (3, 0)$ : $2h_1 + h_2 = 3$ et $h_1 + 2h_2 = 0$, d’où $h_1 = 2$, $h_2 = -1$.

Le nouveau point est $(2, -1)$ : c’est exactement le minimum trouvé au chapitre 5. Comme $f$ est quadratique, son développement de Taylor à l’ordre 2 est exact, et Newton, qui minimise ce développement, tombe sur le minimum en une itération, quel que soit le point de départ.

En code

Python

import numpy as np

grad = lambda v: np.array([2 * v[0] + v[1] - 3, v[0] + 2 * v[1]])
H = np.array([[2.0, 1.0], [1.0, 2.0]])      # Hessienne constante

for depart in ([0.0, 0.0], [10.0, -7.0]):
    x = np.array(depart)
    x = x + np.linalg.solve(H, -grad(x))    # une itération de Newton
    print(depart, "->", x, grad(x))         # toujours [ 2. -1.], gradient nul

R

grad <- function(v) c(2 * v[1] + v[2] - 3, v[1] + 2 * v[2])
H <- rbind(c(2, 1), c(1, 2))                # Hessienne constante

for (depart in list(c(0, 0), c(10, -7))) {
  x <- depart + solve(H, -grad(depart))     # une itération de Newton
  cat(depart, "->", x, "; gradient", grad(x), "\n")
}
# 0 0 -> 2 -1 ; gradient 0 0
# 10 -7 -> 2 -1 ; gradient 0 0 : le minimum en un pas, quel que soit le départ

Exercice · Une direction qui monte

difficile

Pour l’altitude $g(x, y) = x^3 - 3x + y^2$, au point $(-0{,}9 ; 0)$ :

1. Calculez la direction de Newton $h$ et le produit scalaire $\nabla g^T h$. Est-ce une direction de descente ?
2. Même question pour la direction du gradient.
3. Proposez une modification de la Hessienne qui rende la direction de Newton descendante.

Voir la solution

1. $\nabla g(-0{,}9 ; 0) = (3 \times 0{,}81 - 3 ; 0) = (-0{,}57 ; 0)$ et $\nabla^2 g = \mathrm{diag}(-5{,}4 ; 2)$. La direction de Newton vaut $h = (-(-0{,}57)/(-5{,}4) ; 0) \approx (-0{,}1056 ; 0)$. Alors $\nabla g^T h \approx (-0{,}57)(-0{,}1056) \approx +0{,}060 \gt 0$ : c’est une direction de montée. Elle mène vers le col en $(-1, 0)$.
2. La direction du gradient est $-\nabla g = (0{,}57 ; 0)$, et $\nabla g^T(-\nabla g) = -0{,}325 \lt 0$ : elle descend (vers l’est, en direction de la vallée en $(1, 0)$).
3. On remplace la Hessienne par $\nabla^2 g + \lambda I = \mathrm{diag}(-5{,}4 + \lambda ; 2 + \lambda)$, avec $\lambda \gt 5{,}4$ pour la rendre définie positive. Par exemple, avec $\lambda = 6{,}4$, on obtient $\mathrm{diag}(1 ; 8{,}4)$ et la direction $h = (0{,}57 ; 0)$, qui descend. Une matrice définie positive $B$ donne toujours une direction de descente, car $\nabla g^T(-B^{-1}\nabla g) \lt 0$ dès que le gradient est non nul.

En code

Python

import numpy as np

x = np.array([-0.9, 0.0])
g = np.array([3 * x[0]**2 - 3, 2 * x[1]])          # (-0.57, 0)
H = np.diag([6 * x[0], 2.0])                       # diag(-5.4, 2)

h_newton = np.linalg.solve(H, -g)
print(h_newton, g @ h_newton)                      # [-0.1056 0.] 0.0602 > 0 : montée
print(-g, g @ -g)                                  # [0.57 0.] -0.3249 < 0 : descente

h_corrige = np.linalg.solve(H + 6.4 * np.eye(2), -g)   # λ = 6,4 > 5,4
print(h_corrige, g @ h_corrige)                    # [0.57 0.] -0.3249 : descente

R

x <- c(-0.9, 0)
g <- c(3 * x[1]^2 - 3, 2 * x[2])                   # (-0.57, 0)
H <- diag(c(6 * x[1], 2))                          # diag(-5.4, 2)

h_newton <- solve(H, -g)
h_newton; sum(g * h_newton)                        # -0.1056 0 ; 0.0602 > 0 : montée
sum(g * -g)                                        # -0.3249 < 0 : le gradient descend

h_corrige <- solve(H + 6.4 * diag(2), -g)          # λ = 6,4 > 5,4
h_corrige; sum(g * h_corrige)                      # 0.57 0 ; -0.3249 : descente