Optimisation numérique et science des données
Optimisation numérique et science des données
Emmanuel Trélat1
1. Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions (LJLL),
F-75005 Paris, France ([Link]@[Link]).
Table des matières
1 Calcul différentiel 5
1.1 Différentielle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Dérivées d’ordre supérieur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Théorème d’inversion locale, des fonctions implicites . . . . . . . . . . . . . . . . . 10
2 Convexité 12
2.1 Définition et propriétés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Théorème de projection sur un convexe fermé . . . . . . . . . . . . . . . . . . . . . 18
2.3 Sous-différentiabilité des fonctions convexes . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Conjuguée convexe (transformée de Fenchel) . . . . . . . . . . . . . . . . . . . . . . 23
1
4 Minimisation sous contraintes 57
4.1 Existence et unicité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 Conditions d’optimalité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2.1 Conditions d’optimalité du premier ordre sur un ensemble convexe . . . . . 58
4.2.2 Conditions d’optimalité du premier ordre : multiplicateurs de Lagrange . . 58
[Link] Contraintes d’égalité . . . . . . . . . . . . . . . . . . . . . . . . . . 58
[Link] Contraintes d’égalité et d’inégalité : conditions KKT . . . . . . . . 62
[Link] Application : fonctionnelle quadratique avec contraintes d’égalité
affines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2.3 Conditions d’optimalité du second ordre . . . . . . . . . . . . . . . . . . . . 71
[Link] Conditions générales . . . . . . . . . . . . . . . . . . . . . . . . . . 71
[Link] Application à l’analyse de sensibilité . . . . . . . . . . . . . . . . . 73
4.3 Algorithmes d’optimisation avec contraintes . . . . . . . . . . . . . . . . . . . . . . 77
4.3.1 Méthodes primales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
[Link] Méthodes de projection . . . . . . . . . . . . . . . . . . . . . . . . 77
[Link] Méthodes de pénalisation . . . . . . . . . . . . . . . . . . . . . . . 80
[Link] Méthode de Lagrange-Newton . . . . . . . . . . . . . . . . . . . . 84
[Link] Méthode SQP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
[Link] Méthode du Lagrangien augmenté . . . . . . . . . . . . . . . . . . 87
4.3.2 Méthodes duales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
[Link] Point selle du Lagrangien . . . . . . . . . . . . . . . . . . . . . . . 89
[Link] Problème primal et problème dual . . . . . . . . . . . . . . . . . . 91
[Link] Théorème de dualité . . . . . . . . . . . . . . . . . . . . . . . . . . 93
[Link] Méthodes duales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5 Conclusion et compléments 95
5.1 Utilisation de AMPL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.2 Gradient stochastique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.3 Apprentissage, deep learning, rétropropagation . . . . . . . . . . . . . . . . . . . . 104
2
Résumé. Ce cours permet d’acquérir les outils mathématiques théoriques et pratiques de pointe
en optimisation numérique et science des données. L’objectif est d’apprendre à modéliser et résoudre
des problèmes complexes d’optimisation, avec ou sans contraintes, et d’apprendre à mettre en
oeuvre divers algorithmes innovants efficaces pour l’approximation numérique des solutions. Dans
ce cours, on apprendra les méthodes classiques d’optimisation : existence, conditions de premier et
de second ordre, diverses variantes de méthodes de gradient, conditions de Karush-Kuhn-Tucker,
dualité Lagrangienne, puis on fera une ouverture à la science des données : gradient stochastique,
gradient coordonnée par coordonnée, gradient non lisse, TensorFlow. Des TD et TP (en Python)
viendront compléter la formation, ainsi qu’une introduction aux méthodes les plus à la pointe :
différentiation automatique (AMPL) couplée aux outils experts (IpOpt). Elles seront illustrées sur
divers exemples, comme l’analyse d’images ou le machine learning.
Motivation. Les algorithmes de machine learning et deep learning sont basés sur des techniques
d’optimisation. Les problèmes sous-jacents sont souvent posés en très grande dimension, ce qui
rend leur résolution numérique difficile. En grande dimension, les approximations convexes sont
particulièrement utiles, car le problème associé a une solution unique et les algorithmes sont plus
efficaces. Mais en toute généralité, les problèmes sont non convexes et de plus comportent des
contraintes d’égalité ou d’inégalité.
Dans les techniques d’apprentissage, les problèmes de type "moindres carrés" sont d’une grande
importance. Ce sont des problèmes de minimisation convexe
min kAx − bk22
x∈IRn
où A ∈ IRp×n , b ∈ IRp , qui sont souvent posés avec n et/ou p grands. On rencontre beaucoup les
variantes
minn kAx − bk22 + αkxk22
x∈IR
où on pénalise avec la norme `1 , cette dernière ayant l’avantage de promouvoir les solutions par-
cimonieuses (“sparse" en anglais), i.e., la solution x n’a qu’un nombre minimal de composantes
actives. Cette dernière est beaucoup exploitée dans le domaine de l’analyse d’image, depuis les
années 2000 où on a découvert son efficacité et su l’exploiter par des algorithmes d’analyse convexe
et non lisse. Des variantes plus élaborées des moindres carrés comportent aussi des contraintes.
De manière générale, les problèmes de minimisation étudiés dans ce cours s’écrivent
où f : E → IR est une fonction sur un espace vectoriel E topologique (on la prendra souvent de
classe C 1 ou C 2 par la suite, avec E Banach ou Hilbert, mais la fonction pourrait être peu régulière)
et C ⊂ E est un sous-ensemble de E représentant des contraintes. On parle de minimisation sans
contrainte si C = E, et de minimisation sous contraintes lorsque C ( E.
On dit que x∗ ∈ C est un minimiseur global de f sur C si
f (x∗ ) = min f (x).
x∈C
On note aussi
x∗ = argminf (x).
x∈C
3
On dit que x∗ est un minimiseur local de f s’il existe un voisinage ouvert U de x∗ dans E tel que
où
0 if x ∈ C
χC (x) =
+∞ sinon
Cette remarque, formelle pour l’instant, conduit à la notions de pénalisation et aux importantes
méthodes de dualité en optimisation.
4
Chapitre 1
Calcul différentiel
Ce chapitre contient des rappels en calcul différentiel. Le concept essentiel est la notion de
différentielle (de Fréchet, de Gateaux) d’une fonction.
Bien que, en vue de développer des algorithmes de calcul, on s’intéresse prioritairement à des
fonctions f : IRn → IR (donc, sur un espace de dimension finie, éventuellement grande), autant que
possible on s’attachera à donner des énoncés généraux, tant que cela ne complique pas les énoncés.
Ici, tous les espaces vectoriels qu’on va considérer sont réels. On rappelle qu’une norme k k sur
un espace vectoriel E est une application de E dans [0, +∞) qui est homogène (kλxk = |λ|kxk
pour tout λ ∈ IR et tout x ∈ E), sous-additive (inégalité triangulaire kx + yk 6 kxk + kyk pour
tous x, y ∈ E) et qui a la propriété de séparation (pour tout x ∈ E, kxk = 0 ⇒ x = 0).
La notion la plus classique de différentielle, déjà vue dans des cours antérieurs, est la différentielle
de Fréchet d’une fonction f : E → F . Elle s’exprime dans des espaces vectoriels normés. On rappelle
qu’un espace vectoriel normé (E, k kE ) est un Banach s’il est complet. On mentionnera aussi la
notion moins forte de différentielle de Gateaux.
Certaines notions nécessiteront d’avoir un produit scalaire. Or, une norme n’est pas nécessaire-
ment associée à un produit scalaire. Par exemple, dans IRn , la norme k k2 est la norme euclidienne,
issue du produit scalaire euclidien. Mais les normes k k1 et k k∞ ne sont pas des normes associées
à un produit scalaire.
Un espace de Hilbert
p H est un espace de Banach dont la norme découle d’un produit scalaire, par
la formule kxk = hx, xi. On rappelle qu’un produit scalaire est une forme bilinéaire symétrique
et définie positive. Pour tous x, y ∈ H, on a
• kx + yk2 = kxk2 + 2hx, yi + kyk2 (théorème de Pythagore) ;
• |hx, yi| 6 kxkkyk (inégalité de Cauchy-Schwarz).
1.1 Différentielle
Différentielle de Fréchet. Soient E et F des espaces vectoriels normés (en fait, il suffirait que
F soit un espace vectoriel topologique séparé). Soit U un ouvert de E, soit x ∈ U , et soit f : U → F
une application.
On dit que f est différentiable au sens de Fréchet en x s’il existe une application linéaire continue
df (x) : E → F telle que, pour tout h ∈ E tel que x + h ∈ U ,
5
Autrement dit, f est Fréchet-différentiable en x si et seulement si f a un développement limité à
l’ordre 1 en x.
L’application linéaire continue df (x) ∈ L(E, F ) s’appelle la différentielle de Fréchet en x.
Notons que si f est Fréchet différentiable en x alors elle est continue en x.
f (x + th) − f (x)
f 0 (x; h) = lim
t→0
t6=0
t
existe. On parle de dérivée directionnelle au sens de Dini si, dans la limite ci-dessus, on impose de
plus t > 0. On dit que f est différentiable au sens de Gateaux en x si sa dérivée directionnelle en
x existe selon toute direction et si l’application h 7→ f 0 (x; h) est linéaire continue.
Avec un abus de notation, on note aussi df (x).h = f 0 (x; h).
a des dérivées directionnelles en (0, 0) dans toutes les directions, mais elle n’est pas Gateaux
différentiable car l’application h 7→ f 0 ((0, 0), h) n’est pas linéaire.
Dans la suite, quand on parlera de différentielle sans préciser, cela sous-entend différentielle au
sens de Fréchet.
Remarque 2. Pour f : IR → IR, la dérivée usuelle est définie par
f (x + h) − f (x)
f 0 (x) = lim
h→0
h6=0
h
6
• La différentielle au sens de Fréchet en a ∈ IRn \ {0} de la fonction f : IRn → IR définie
a
par f (x) = kxk2 est df (a).h = h kak 2
, hi (produit scalaire euclidien). La fonction f n’est pas
différentiable en 0.
• Soit A ∈ Mn (IR) une matrice symétrique. La différentielle au sens de Fréchet en a ∈ IRn de
l’application f : IRn → IR définie par f (x) = x> Ax est df (a).h = a> Ah + h> Aa = 2a> Ah.
Dérivées partielles. Comme dit ci-dessus, si f est différentiable au sens de Fréchet en x (ou,
seulement, Gateaux différentiable en x) et si h est un vecteur de E, la dérivée directionnelle de f
dans la direction h est
f (x + th) − f (x)
f 0 (x; h) = df (x).h = lim .
t→0 t
Lorsque E = IRn et F = IR, en notant (e1 , . . . , en ) la base canonique, on a
∂f f (x1 , . . . , xj−1 , xj + t, xj+1 , . . . , xn ) − f (x)
(x) = df (x).ej = f 0 (x; ej ) = lim
∂xj t→0
t6=0
t
Notons que l’existence de dérivées partielles n’implique pas forcément la différentiabilité. Par
exemple la fonction f : IR2 → IR définie par f (x, y) = x2xy
+y 2 si (x, y) 6= (0, 0) et f (0, 0) = 0 admet
des dérivées partielles nulles en (0, 0), mais n’est pas différentiable en (0, 0) car elle n’y est pas
continue.
On a df (x).h = Jf (x)h.
Un cas particulier important est lorsque p = 1, i.e., pour une fonction f : IRn → IR. En notant
h , i le produit scalaire euclidien de IRn , on a
où ∂f
∂x1 (x)
∇f (x) =
..
.
∂f
∂xn (x)
est le gradient de f en x.
De manière plus générale, le gradient est défini dès que E est un espace de Hilbert (il faut un
produit scalaire) : pour une fonction f : E :→ IR différentiable, avec E Hilbert, pour tout x ∈ E
le vecteur ∇f (x) est l’unique vecteur de E tel que df (x).h = h∇f (x), hi pour tout h ∈ E, où h , i
est le produit scalaire sur E.
7
Propriétés.
• Si f et g sont deux applications différentiables en x, alors pour tous λ, µ ∈ IR, l’application
λf + µg est aussi différentiable en x, et
df : U → L(E, F )
x 7→ df (x)
On suppose que E et F sont des espaces vectoriels normés, de sorte que l’espace L(E, F ) des
applications linéaires continues de E dans F est un espace vectoriel normé, équipé de la norme
d’opérateur définie par
k`(x)kF
k`kL(E,F ) = sup = sup k`(x)kF ∀` ∈ L(E, F ).
x∈E\{0} kxkE x∈E
kxkE =1
On dit que f est de classe C 1 en x si elle est différentiable sur un voisinage ouvert U de x et si
l’application df : U → L(E, F ) est continue en x.
8
Théorème 2. (Théorème de Schwarz) Si f est deux fois différentiable en x alors d2 f (x) est une
application bilinéaire symétrique, i.e.,
Un cas particulier important est lorsque E = IRn et F = IR. Alors d2 f (x) est une forme
bilinéaire symétrique (donc, associée à une forme quadratique) qu’on représente par une matrice
appelée la Hessienne de f en x, qui est la matrice symétrique
∂2f
Hf (x) = (x)
∂xi ∂xj 16i,j6n
et on a
d2 f (x).(h, h) = h> Hf (x)h ∀h ∈ IRn
De façon générale, par itération, on dit que f : E → F est k fois différentiable en x si elle
est (k − 1) fois différentiable sur un voisinage ouvert U de a et si sa différentielle (k − 1)ème
dk−1 f : U → L(E k−1 , F ) est différentiable en x. On note dk f (x).(h1 , . . . , hk ) l’action de dk f (x)
sur (h1 , . . . , hk ) ∈ E k . Le théorème de Schwarz se généralise à l’ordre k : si f est k fois différentiable
en x alors dk f (x) est une application k-linéaire symétrique.
On dit que f est de classe C k si l’application dk f est continue.
Formules de Taylor. On rappelle d’abord la formule de Taylor avec reste intégral. Soit U est
un ouvert de E, soit f : U → F une application de classe C k+1 , soient x ∈ U et h ∈ E tels que le
segment [x, x + h] est contenu dans U . Alors
1 1
f (x + h) = f (x) + df (x).h + d2 f (x).(h, h) + · · · + dk f (x).(h, . . . , h)
2 k!
1 1
Z
+ (1 − t)k dk+1 f (x + th).(h, . . . , h) dt
k! 0
Ce résultat est très facile à obtenir en se ramenant à la dimension 1, le long du segment [x, x + h].
Comme nous allons beaucoup utiliser cette technique par la suite, nous la rappelons ici. On pose
ϕ(t) = f (x + th)
R1
de façon à ce que ϕ(0) = f (x) et ϕ(1) = f (x+h). On part de ϕ(1) = ϕ(0)+ 0 ϕ0 (t) dt. En intégrant
R1
par parties, on dérive ϕ0 et on intègre 1 en −(1−t), on obtient ϕ(1) = ϕ(0)+ϕ0 (0)+ 0 (1−t)ϕ00 (t) dt,
puis en itérant,
1 1
Z
1
ϕ(1) = ϕ(0) + ϕ0 (0) + · · · + ϕ(k) (0) + (1 − t)k ϕ(k+1) (t) dt
k! k! 0
Il reste à noter que ϕ0 (0) = df (x).h, . . . , ϕ(k) (0) = dk f (x).(h, . . . , h) et ϕ(k+1) (t) = dk+1 f (x +
th).(h, . . . , h).
Rappelons la formule de Taylor avec reste en o. Si f est k fois différentiable en x alors
1 1
f (x + h) = f (x) + df (x).h + d2 f (x).(h, h) + · · · + dk f (x).(h, . . . , h) + o(khkkE )
2 k!
9
lorsque h → 0. Attention, pour k > 1, f peut avoir un développement limité d’ordre k en x, et
pour autant, ne pas être k fois différentiable en x. L’équivalence n’est vraie que pour k = 1.
Un cas particulier important est lorsque E = IRn et k = 2 : si f : IRn → IR est deux fois
différentiable en x alors
1
f (x + h) = f (x) + df (x).h + d2 f (x).(h, h) + o(khk2 )
2
1
= f (x) + ∇f (x) h + h> Hf (x)h + o(khk2 )
>
2
Exemple 2. Un exemple qu’on utilisera par la suite est la fonction f : IRn → IR définie par
1
f (x) = x> Ax − b> x ∀x ∈ IRn
2
où A est une matrice symétrique réelle, et b ∈ IRn . On a ici
df (x).h = x> Ah − b> h, i.e., ∇f (x) = Ax − b et Hf (x) = A
10
Théorème des fonctions implicites. L’objectif est de résoudre une équation f (x, y) = 0, où
f : E × F → G, et d’exprimer y comme fonction de x (localement). Le théorème des fonctions
implicites est en fait équivalent au théorème d’inversion locale.
Par exemple, en dimension finie, supposons que f soit linéaire : f (x, y) = A1 x + A2 y où A1 et
A2 sont des matrices. Pour résoudre A1 x + A2 y = 0 en exprimant y comme fonction de x, on a
besoin de supposer que A2 est inversible, et alors on trouve y = −A−1 2 A1 x. Dans le cas général
non linéaire, A2 est la différentielle de f par rapport à y et on a le théorème suivant.
Théorème 4. Soient E, F et G des Banach. Soit f : E×F → G de classe C 1 , et soit (a, b) ∈ E×F
tel que f (a, b) = 0. On suppose que ∂f∂y (a, b) ∈ L(F, G) (qui est la différentielle de f , uniquement
par rapport à la variable y ∈ F ) est inversible. Alors il existe un voisinage ouvert V de a dans E,
un voisinage ouvert W de b dans F , et une (unique) application ϕ : V → W de classe C 1 telle que
f (x, ϕ(x)) = 0 pour tout x ∈ V , avec ϕ(a) = b.
Autrement dit, on a résolu localement autour de (a, b) l’équation f (x, y) = 0 par rapport à x
(en exprimant, localement, y comme fonction de x).
Si f est de classe C k avec k > 1 alors ϕ est C k .
Pour exprimer la différentielle de ϕ, il suffit de dériver par rapport à x l’égalité f (x, ϕ(x)) = 0.
On obtient ∂f ∂f
∂x (x, ϕ(x)) + ∂y (x, ϕ(x)) ◦ dϕ(x) = 0, d’où
−1
∂f ∂f
dϕ(x) = − (x, ϕ(x)) (x, ϕ(x)).
∂y ∂x
f (a + δx, b + δy) = 0
∂f ∂f
(a, b).δx + (a, b).δy ' 0
∂x ∂y
qui est exactement δy = dϕ(a).δx, comme attendu puisque, en toute rigueur on sait que b + δy =
ϕ(a + δx).
11
Chapitre 2
Convexité
Autrement dit, le segment [x, y] est inclus dans C, pour tous points x, y ∈ C.
Dans IR, les ensembles convexes sont exactement les intervalles. Une union disjointe d’intervalles
n’est pas convexe.
Définition 2. Soit C ⊂ E un sous-ensemble convexe de l’espace vectoriel E. Une fonction f :
C → IR ∪ {+∞} est dite convexe si
autrement dit la valeur prise par f sur toute combinaison convexe de points xi est plus petite que
la combinaison convexe des valeurs de f aux points xi .
12
On dit que f est strictement convexe si
On suppose que E est un espace vectoriel normé. Pour α > 0, la fonction f est dite α-convexe
(ou fortement convexe de module α) si
α
∀x, y ∈ C ∩ Dom(f ) ∀λ ∈ [0, 1] f (λx + (1 − λ)y) 6 λf (x) + (1 − λ)f (y) − λ(1 − λ)kx − yk2E
2
Dans la littérature existante, on trouve parfois le vocabulaire suivant à propos des fonctions
convexes : lorsque Dom(f ) 6= ∅, on dit que f est propre.
On a les implications : α-convexe ⇒ strictement convexe ⇒ convexe.
Exemple 3. • Dans un espace vectoriel normé, la fonction f (x) = kxk est convexe.
n
Dans IR , la fonction f (x) = kxk22 est strictement convexe. La fonction f (x) = kxk2∞ ou kxk21
est convexe mais pas strictement convexe.
• Toute application linéaire ou affine est convexe (mais pas strictement convexe).
• Soit A ∈ Mn (IR) une matrice carrée symétrique positive d’ordre n. Soit b ∈ IRn . La fonction
1 >
f (x) = x Ax − b> x
2
est convexe. Elle est strictement convexe lorsque A est définie positive.
Plus généralement, dans un Hilbert H, soit A : H → H un opérateur linéaire autoadjoint.
On suppose que A est monotone, i.e., hA(x − y), x − yiH > 0 pour tous x, y ∈ H (de manière
équivalente, puisque A est linéaire : hAx, xi > 0 pour tout x ∈ H ; on dit aussi que −A est
dissipatif). Soit b ∈ H. La fonction
1
f (x) = hx, AxiH − hb, xiH
2
est convexe.
13
Remarque 7. Une fonction f est convexe si et seulement si son épigraphe
Si une fonction est convexe, alors, pour tout α ∈ IR, l’ensemble de “sous-niveau"
{x ∈ Dom(f ) | f (x) 6 α}
(“sublevel
p set" en anglais) est convexe. La réciproque est fausse : par exemple la fonction f (x) =
|x| n’est pas convexe, mais ses ensembles de sous-niveau sont convexes (on dit qu’elle est “quasi-
convexe").
Remarque 8. Pour les fonctions f : IR → IR, on rappelle les propriétés suivantes :
• lorsque f est dérivable, f est (strictement) convexe si et seulement si f 0 est (strictement)
croissante, si et seulement si le graphe de f est (strictement, sauf au point de tangence)
au-dessus de ses tangentes ;
• lorsque f est deux fois dérivable, f est convexe si et seulement si f 00 > 0 ; f est strictement
convexe si et seulement si f 00 > 0 et ne s’annule que sur un ensemble d’intérieur vide.
(NB : f (x) = x4 est strictement convexe mais f 00 (0) = 0)
Remarque 9. Comme on l’a déjà vu pour la formule de Taylor, on peut démontrer de nombreuses
propriétés en multi-D à partir des propriétés en 1D. La technique est toujours la même : étant
donnée une fonction f : E → IR (où E est un espace vectoriel), étant donnés deux points x, y ∈ E,
on regarde la fonction f le long du segment [x, y] en posant
14
Continuité des fonctions convexes.
Théorème 5. Soit E un espace vectoriel topologique, et soit f : E → IR ∪ {+∞} une fonction
convexe. Si f est localement bornée autour d’un point x ∈ E (qui est dans l’intérieur de Dom(f ))
alors f est continue en x.
Comme conséquence, si E est de dimension finie alors f est continue sur l’intérieur de Dom(f ).
Démonstration. Sans perte de généralité, on suppose que 0 est dans l’intérieur de Dom(f ) et que
f (0) = 0. Soit V un voisinage ouvert convexe de 0 et soit a > 0 tel que f (x) 6 a pour tout x ∈ V .
L’ouvert W = V ∩ (−V ) est aussi un voisinage ouvert convexe de 0, qui est de plus symétrique par
rapport à 0). Soit ε ∈]0, 1[.
Pour tout x ∈ εW (i.e., xε ∈ W ), x s’écrit comme la combinaison convexe x = (1 − ε)0 + ε xε ,
donc par convexité de f ,
x x
f (x) 6 (1 − ε)f (0) + εf = εf 6 εa.
ε ε
ε −x
Par ailleurs, on a − xε ∈ W , et comme 0 s’écrit comme la combinaison convexe 0 = 1+ε
1
x + 1+ε ε ,
par convexité de f ,
1 ε −x
0 = f (0) 6 f (x) + f
1+ε 1+ε ε
−x
donc f (x) > −εf ε > −εa. De ces deux inégalités, on déduit que |f (x)| 6 εa pour tout x ∈ εW .
On conclut que f est continue en 0.
En dimension finie, si l’intérieur de Dom(f ) est non vide, il contient n + 1 points xi affinement
indépendants. Par convexité, on trouve que
n+1 n+1
! n+1
X X X
∀λ1 , . . . , λn+1 > 0 | λi = 1 f λi xi 6 λi f (xi ) 6 max f (xi ).
16i6n+1
i=1 i=1 i=1
15
On parle de minorantes affines de f .
16
Caractérisation de la convexité par la Hessienne.
Théorème 7. Soit E un espace vectoriel normé et soit f : E → IR une fonction deux fois
différentiable sur Dom(f ). La fonction f est convexe si et seulement si, pour tout x ∈ Dom(f ),
d2 f (x) est une forme quadratique positive. Si d2 f (x) est définie positive en tout x ∈ Dom(f ) alors
f est strictement convexe (la réciproque est fausse : prendre f (x) = x4 ).
Lorsque E est de dimension finie, cela s’exprime sous la forme Hf (x) > 0 , i.e., f est convexe
si et seulement si sa Hessienne est en tout point une matrice symétrique positive. Si la Hessienne
est définie positive en tout point alors f est strictement convexe.
Démonstration. Suivant la remarque 9, en posant ϕ(t) = f (x + t(x − y)), f est convexe si et
seulement si ϕ est convexe pour tous x, y ∈ E. On a ϕ00 (t) = d2 f (x + t(y − x)).(y − x, y − x). Or,
en 1D, ϕ est convexe si et seulement si ϕ00 > 0. Le résultat s’ensuit.
Exemple 4. Soit A ∈ Mn (IR) une matrice symétrique réelle et soit b ∈ IRn . La fonction
1 >
f (x) = x Ax − b> x
2
est convexe (resp., strictement convexe) si et seulement si A > 0 (resp., A > 0).
Théorème 8. Soit H un Hilbert (espace préhilbertien suffirait), soit C ⊂ H un sous-ensemble
convexe et soit f : C → IR une fonction différentiable sur C. Soit α > 0.
f est α-convexe sur C
⇔ x 7→ f (x) − α2 kxk2 est convexe sur C
⇔ ∀x, y ∈ C f (y) > f (x) + df (x).(y − x) + α2 ky − xk2
⇔ ∀x, y ∈ C (df (x) − df (y)).(x − y) > αkx − yk2
Si f est deux fois différentiable, alors f est α-convexe si et seulement si d2 f (x).(h, h) > αkhk2 pour
tout x ∈ C, h ∈ H. En dimension finie, cette condition s’exprime sous la forme Hf (x) > αIn .
La troisième propriété stipule que non seulement le graphe de f est au-dessus de ses tangentes,
mais que, en plus, entre les deux on peut placer une parabole. Cela veut donc dire que le graphe
de f s’éloigne (au-dessus) de ses tangentes au moins comme un carré.
Démonstration. Posons g(x) = f (x) − α2 kxk2 . La fonction g est convexe si et seulement si tg(x) +
(1 − t)g(y) − g(tx + (1 − t)y) > 0 pour tout t ∈ [0, 1] et tous x, y ∈ C. Or,
et ainsi on trouve que g est convexe si et seulement si f est α-convexe. Le reste est facile.
17
2.2 Théorème de projection sur un convexe fermé
Soit H un espace de Hilbert. On rappelle tout d’abord l’identité du parallélogramme :
2 2
x+y x−y 1
kxk2 + kyk2
∀x, y ∈ H + =
2 2 2
donc la suite est de Cauchy, donc elle converge (car l’espace est complet) vers un x∗ ∈ C (car C
est fermé), et on obtient (2.2).
Montrons que (2.2) ⇒ (2.1). Pour tout t ∈]0, 1], on a, par la propriété de minimisation,
18
donc en développant, 0 6 2thx − x∗ , x∗ − zi + t2 kx∗ − zk2 , puis on divise par t et on fait tendre
t → 0 et on obtient (2.1).
Montrons que (2.1) ⇒ (2.2). On a
d’où (2.2).
Montrons que PC est 1-Lipschitzienne. On note x∗ = PC (x) et y ∗ = PC (y). Par (2.1), on a
hx − x∗ , z − x∗ i 6 0 et hy − y ∗ , z − y ∗ i 6 0 ∀z ∈ C.
kx∗ − y ∗ k2 6 hx − y, x∗ − y ∗ i ∀x, y ∈ H
19
2.3 Sous-différentiabilité des fonctions convexes
Soit H un espace de Hilbert.
Définition 3. Soit f : H → IR ∪ {+∞} une fonction convexe. Un vecteur p ∈ H est appelé
sous-gradient de f au point x ∈ Dom(f ) si
Lorsque la fonction convexe f est différentiable au sens de Fréchet en x ∈ H, alors ∂f (x) est le
singleton :
∂f (x) = {∇f (x)}.
Lemme 1. Pour tout x ∈ Dom(f ), ∂f (x) est un ensemble convexe fermé.
Démonstration. Soient p1 , p2 ∈ ∂f (x). Alors, pour tout y ∈ Dom(f ), f (y) > f (x) + hp1 , y − xi et
f (y) > f (x) + hp2 , y − xi, d’où, pour tout λ ∈ [0, 1], f (y) > f (x) + hλp1 + (1 − λ)p2 , y − xi, et donc
λp1 + (1 − λ)p2 ∈ ∂f (x). Le caractère fermé est évident.
La notion de sous-différentiel n’est pas restreinte aux espaces de Hilbert. Lorsque H est un
espace vectoriel topologique, on remplace dans la définition le produit scalaire par le crochet de
dualité h , iH 0 ,H .
La notion de sous-différentiel n’est pas restreinte aux fonctions convexes. On peut définir la
sous-différentiabilité d’une fonction dans un cadre très général (théorie d’analyse non lisse). Mais,
dans cette courte section introductive, on se limite aux fonctions convexes.
20
Définition 4. Pour E espace topologique, une fonction f : E → IR est semi-continue inférieure-
ment (sci) en x ∈ E si pour tout ε > 0 il existe un voisinage ouvert U de x dans E tel que, pour
tout y ∈ U , f (x) 6 f (y) + ε, ou de manière équivalente,
La fonction f est sci sur E si et seulement si, pour tout α ∈ IR, l’ensemble {x ∈ E | f (x) 6 α}
est fermé, si et seulement si l’épigraphe Epi(f ) = {(x, α) ∈ Dom(f ) × IR | f (x) 6 α} est fermé.
Bien que ce ne soit pas utile ici, précisons qu’on peut définir, de même, la notion de fonction
semi-continue supérieurement (scs) : simplement, f est scs si et seulement si −f est sci (i.e.,
f (x) > lim supy→x f (y)).
Théorème 10. Toute fonction convexe sci f : H → IR ∪ {+∞} est sous-différentiable sur l’inté-
rieur de son domaine : pour tout x appartenant à l’intérieur de Dom(f ), l’ensemble ∂f (x) est un
convexe non vide. Il est de plus borné si H est de dimension finie.
Démonstration. Comme f est convexe sci, Epi(f ) est un ensemble convexe fermé. Pour x ∈
Dom(f ), le point (x, f (x)) appartient à la frontière de Epi(f ). Donc par le théorème de sépa-
ration convexe (voir remarque 11 et voir la figure), il existe (p, p0 ) ∈ H × IR \ {(0, 0)} (défini à
scalaire multiplicatif > 0 près) tel que
p y−x
, 60 ∀(y, α) ∈ Epi(f )
p0 α − f (x)
de sorte que
hp, y − xi + p0 (α − f (x)) 6 0.
En prenant y = x et α = f (x) + t, t > 0, on voit que, forcément, p0 6 0.
α
y−x
α − f (x)
f (x) p
p0
x y
Montrons par l’absurde que p0 < 0. Si p0 = 0 alors hp, y − xi 6 0 pour tout y ∈ Dom(f ).
Comme x est dans l’intérieur de Dom(f ), il existe ε > 0 tel que la boule fermée B̄(x, ε) soit incluse
dans l’intérieur de Dom(f ). En prenant y = x + h avec h ∈ B̄(x, ε), on obtient donc hp, hi = 0
pour tout h ∈ B̄(x, ε). Mais cela implique p = 0 (en effet, prendre h = tp avec t positif et négatif
et |t| assez petit). Autrement dit, (p, p0 ) = (0, 0), ce qui est une contradiction car (p, p0 ) 6= (0, 0).
Comme le couple (p, p0 ) peut être multiplié par un scalaire strictement positif sans changer
les inégalités ci-dessus, quitte à le multiplier par −1/p0 on se ramène à p0 = −1. Ainsi, on a
21
obtenu hp, y − xi 6 α − f (x) pour tout (y, α) ∈ Epi(f ). En particulier, pour α = f (y), cela donne
f (y) > f (x) + hp, y − xi. On a bien obtenu l’existence d’un p ∈ ∂f (x), i.e., ∂f (x) 6= ∅.
Il reste à montrer que ∂f (x) est borné si H est de dimension finie. Par l’absurde, s’il n’est pas
borné alors il existe une suite (pk )k∈IN d’éléments de ∂f (x), tels que kpk k → +∞. Par définition,
on a f (y) > f (x) + hpk , y − xi pour tout y ∈ Dom(f ), pour tout k ∈ IN. Divisons par kpk k et
posons Ψk = kppkk k . On a
1 1
f (y) > f (x) + hΨk , y − xi ∀y ∈ Dom(f ).
kpk k kpk k
Le vecteur Ψk appartient à la sphère unité de H qui est compacte car H est de dimension finie,
donc à sous-suite près on a Ψk → Ψ avec kΨk = 1. Comme kpk k → +∞, on obtient en passant à
la limite hΨ, y − xi = 0 pour tout y ∈ Dom(f ). Comme x est dans l’intérieur de Dom(f ), on en
déduit comme précédemment que Ψ = 0, ce qui est absurde puisque kΨk = 1.
La notion de sous-différentiabilité généralise la notion de différentiabilité (y compris pour des
fonctions non convexes). Lorsqu’une fonction f : H → IR est différentiable au sens de Fréchet en
x ∈ H, alors ∂f (x) = {∇f (x)}. Il existe une théorie très développée de la sous-différentiabilité
(appelée “analyse non lisse"), et beaucoup de propriétés des fonctions sous-différentiables, avec tout
un “calcul sous-différentiel" (modules d’analyse convexe ou analyse non lisse, généralement étudiés
au niveau M2).
Ici, on se contente d’une brève introduction, en se restreignant de plus aux fonctions convexes,
mais on mentionne le fait très simple suivant : le sous-différentiel permet notamment de caractériser
les minimiseurs des fonctions convexes.
Théorème 11. Soit f : H → IR une fonction convexe sous-différentiable en x∗ ∈ H. Le point
x∗ ∈ H est un minimiseur de la fonction f sur H si et seulement si 0 ∈ ∂f (x∗ ) :
pour tout x ∈ Dom(f ). Réciproquement, si f (x) > f (x∗ ) pour tout x ∈ Dom(f ) alors 0 ∈ ∂f (x∗ )
par définition.
Remarque 13. Soit f : H → IR ∪ {+∞} une fonction convexe telle que Dom(f ) 6= ∅. Soit
x ∈ Dom(f ). Les propriétés suivantes sont équivalentes :
• p ∈ H est un sous-gradient de f en x ;
• ∀y ∈ H f (y) > f (x) + hp, y − xi ;
• ∀h ∈ H f 0 (x; h) > hp, hi ;
• x est un minimiseur de la fonction y ∈ H 7→ f (y) − hp, yi ;
• f (x) + f ∗ (p) > hp, xi ;
• f (x) + f ∗ (p) = hp, xi.
En particulier,
f 0 (x; h) = sup hp, hi
p∈∂f (x)
Dans les deux derniers items, en anticipant légèrement, on a utilisé la conjuguée convexe f ∗ , définie
dans la section suivante.
22
2.4 Conjuguée convexe (transformée de Fenchel)
La fonction conjuguée, appelée aussi transformée de Fenchel, est utilisée pour calculer le sous-
différentiel d’une fonction convexe, caractériser des problèmes duaux (voir plus loin dans ce cours),
ou encore, en utilisant la biconjuguée, pour convexifier une fonction non convexe (voir ci-dessous).
Définition 5. Soit H un espace de Hilbert. Soit f : H → IR ∪ {+∞} une fonction quelconque. La
fonction conjuguée (ou transformée de Fenchel) f ∗ de f est définie par
La conjuguée f ∗ est une fonction sur H à valeurs dans IR ∪ {+∞}. Son domaine est défini par
Dom(f ∗ ) = {p ∈ H | f ∗ (p) < +∞}.
De même que pour le sous-différentiel, la notion de conjuguée convexe se généralise lorsque H
est un espace vectoriel topologique. Dans ce cas, la fonction f ∗ est définie sur H 0 , et le produit
scalaire ci-dessus est remplacé par le crochet de dualité.
Par définition, on a l’inégalité de Fenchel
hp, xi − β 6 f (x) ∀x ∈ H ⇔ β > hp, xi − f (x) ∀x ∈ H ⇔ β > sup (hp, xi − f (x)) = f ∗ (p).
x∈H
23
Autrement dit, à p fixé, l’hyperplan affine x 7→ hp, xi − f ∗ (p) est la plus grande minorante affine
de f . Cet hyperplan a un point de contact avec le graphe de f (cf la proposition 1 ci-dessous).
q0
Exemple 5. Pour f (x) = 1q kxkq , on a f ∗ (p) = 1
q 0 kxk pour q ∈]1, +∞[ et 1
q + 1
q0 = 1. L’inégalité
1 q 1 q0
de Fenchel donne alors hx, pi 6 q kxk
+ q 0 kpk
pour tous x, p ∈ H, ce qui est l’inégalité de
Young généralisée (bien connue pour q = q 0 = 2). Ainsi, l’inégalité de Fenchel permet d’obtenir de
nouvelles inégalités.
La propriété suivante montre que les concepts de sous-différentiel et de transformée de Fenchel
sont étroitement liés, et caractérise le cas d’égalité dans l’inégalité de Fenchel (on n’a pas besoin
de supposer f convexe, à condition de garder toutefois la même définition d’un sous-gradient).
Démonstration. On a
p ∈ ∂f (x) ⇔ ∀y ∈ H f (y) > f (x) + hp, y − xi
⇔ ∀y ∈ H hp, yi − f (y) 6 hp, xi − f (x) avec égalité pour y = x
⇔ sup (hp, yi − f (y)) = hp, xi − f (x)
y∈H
| {z }
f ∗ (p)
Comme f ∗ , la fonction f ∗∗ est une fonction sur H à valeurs dans IR ∪ {+∞}, et elle est convexe
sci (comme supremum de fonctions affines, ou bien, parce que c’est une conjuguée).
Théorème 12. (Fenchel-Moreau) On a toujours
f ∗∗ 6 f
et en fait, la biconjuguée f ∗∗ est la plus grande fonction convexe sci inférieure ou égale à f . On a
24
Démonstration. On définit l’ensemble S des minorantes affines de f :
S = {(p, β) ∈ H × (IR ∪ {+∞}) | ∀x ∈ H f (x) > hp, xi − β}
= {(p, β) ∈ H × (IR ∪ {+∞}) | β > sup (hp, xi − f (x))}
x∈H
∗
= {(p, β) ∈ H × (IR ∪ {+∞}) | β > f (p)}
Pour tout x ∈ H fixé, on a, par définition, hp, xi − β 6 f (x) pour tout (p, β) ∈ S , donc en passant
au sup, on a
sup (hp, xi − β) 6 f (x).
(p,β)∈S
En fait, ce sup vaut exactement f (x), ce qui montre donc que f ∗∗ (x) 6 f (x). En effet :
∗∗
(on a utilisé que le supremum des hp, xi − β sur l’ensemble des β tels que −β 6 −f ∗ (p) est égal à
hp, xi − f ∗ (p)). Pour montrer que f ∗∗ est la plus grande fonction convexe sci inférieure à f , il suffit
de montrer la propriété sur les épigraphes. Pour tout (p, β) ∈ S , comme x 7→ hp, xi − β est une
minorante affine de f , son épigraphe est un demi-espace fermé contenant Epi(f ). Comme, d’après
(2.3), la fonction convexe f ∗∗ est le sup sur (p, β) ∈ S de ces minorantes affines, son épigraphe
Epi(f ∗∗ ) est l’intersection de tous ces demi-espaces fermés contenant Epi(f ), autrement dit, c’est
exactement Conv(Epi(f )) (voir remarque 12).
Epi(f )
f ∗∗
25
Chapitre 3
min f (x)
x∈E
pour tout x ∈ U .
lim f (x) = +∞
kxk→+∞
On dit aussi parfois que f est propre (par définition, une application est propre si l’image
réciproque de tout compact est un compact), mais comme on l’a vu le mot propre peut avoir
plusieurs significations. On préfère ici utiliser le vocabulaire “infinie à l’infini".
Par exemple, f (x) = kxk est infinie à l’infini. La fonction f : IRn → IR définie par f (x) =
1 > >
2 x Ax − b b où A est une matrice symétrique définie positive est infinie à l’infini. Mais par
exemple la fonction f : IR2 → IR définie par f (x) = x21 ne l’est pas.
Théorème 13. Soit f : IRn → IR une fonction continue infinie à l’infini. Alors le problème
min f (x)
x∈IRn
26
admet au moins un minimiseur (global), i.e., il existe au moins un point x∗ ∈ IRn tel que f (x∗ ) =
minx∈IRn f (x). On note aussi x∗ = argminx∈IRn f (x).
Démonstration. Soit A > inf x∈IRn f (x). Comme f est infinie à l’infini, il existe R > 0 tel que
f (x) > A pour tout x ∈ B(0, R) = {x ∈ IRn | kxk > R}. Cela implique que inf x∈IRn f (x) =
inf x∈B(0,R) f (x). Or, comme on est en dimension finie, B(0, R) est compacte, donc f , qui est
continue, atteint son minimum.
Il n’y a pas unicité du minimiseur en général. On peut aussi avoir des minimiseurs locaux qui
ne sont pas globaux. Pour avoir unicité du minimiseur global, un bon moyen est de supposer la
convexité.
Théorème 14. Soit E un espace vectoriel. Si f : E → IR est strictement convexe alors elle a au
plus un minimiseur (forcément global).
Démonstration. En effet, par l’absurde, s’il y avait deux minimiseurs x1 et x2 (i.e., f (x1 ) = f (x2 ) =
minx∈E f (x)), alors, par convexité stricte, on aurait f ( x1 +x
2
2
) < 21 (f (x1 ) + f (x2 )) = minx∈E f (x)
ce qui est absurde.
Par contre, la convexité (même stricte) de f n’implique pas forcément l’existence d’un mini-
miseur. Par exemple, f (x) = ex est strictement convexe mais n’a pas de minimiseur. Son infimum
sur IR vaut 0 mais n’est pas atteint.
Corollaire 2. Si f : IRn → IR est continue, infinie à l’infini et strictement convexe, alors il existe
un unique minimiseur (forcément global) de f .
Remarque 14. Si f : IRn → IR est α-convexe (avec α > 0) alors f est infinie à l’infini (appliquer
le théorème 8), continue et strictement convexe, donc elle admet un unique minimiseur (forcément
global).
En pratique, lorsqu’on minimise une fonction f sur IRn , il est rare que celle-ci soit convexe.
Mais, souvent, on est capable de dire qu’elle est au moins dans une région donnée, et alors, dans
cette région, on arrive à caractériser un unique minimiseur. C’est sur cette vision “locale" que se
basent les algorithmes qu’on verra par la suite et qui sont basés sur des conditions nécessaires
d’optimalité du premier ordre (conditions de dérivée nulle).
Obtenir numériquement un minimiseur global, alors qu’il existe plusieurs (éventuellement, beau-
coup de) minimiseurs locaux, est difficile et fait appel à d’autres techniques (“optimisation globale")
qui sortent du cadre de ce cours.
27
En général, un point qui annule la différentielle de f s’appelle un point extremum. Ainsi, tout
minimiseur local est un point extremum.
Notons bien qu’il s’agit d’une condition nécessaire d’optimalité. Le fait que la dérivée de f
s’annule en x∗ n’implique pas que x∗ soit un minimiseur (local) : ce serait être un maximiseur,
ou bien, un extremum qui n’est ni un minimum ni un maximum (comme l’est 0 pour la fonction
f (x) = x3 ).
La condition est nécessaire et suffisante dans le cas convexe.
Théorème 16. Soit f : E → IR une fonction convexe Gateaux différentiable. Alors x∗ ∈ E est un
minimiseur global si et seulement si df (x∗ ) = 0.
Plus généralement, on a déjà vu que si f est convexe et sous-différentiable, alors x∗ est mini-
miseur global de f si et seulement si 0 ∈ ∂f (x∗ ).
Démonstration. Comme f est convexe, le graphe de f est au-dessus de ses tangentes : f (x∗ + h) >
f (x∗ ) + df (x∗ ).h pour tout h ∈ E, donc si df (x∗ ) = 0 alors x∗ est un minimiseur global de f .
Remarque 15. Même pour f convexe, on n’a pas forcément unicité du minimiseur global : la
fonction f pourrait en effet être constante au voisinage de x∗ .
Lorsque f est convexe, l’ensemble des minimiseurs (globaux) est un sous-ensemble convexe de
E. En effet, si x et y sont deux minimiseurs de f , i.e., f (x) = f (y) = min f , alors pour tout
λ ∈ [0, 1], f (λx + (1 − λ)y) 6 λf (x) + (1 − λ)f (y) = min f donc en fait on a égalité et λx + (1 − λ)y
est aussi un minimiseur.
(on dit que d2 f (x∗ ) est une forme quadratique coercive, ou α-elliptique) alors x∗ est un
minimiseur local strict de f .
Minimiseur local strict signifie que f (x∗ + h) > f (x∗ ) pour tout h 6= 0 de norme assez petite.
28
Montrons maintenant la condition suffisante. On fait un développement de Taylor de f au
second ordre en x∗ , et comme df (x∗ ) = 0, on a
1
f (x∗ + h) = f (x∗ ) + d2 f (x∗ ).(h, h) + o(khk2 )
2
lorsque h → 0. Comme d2 f (x∗ ) est α-elliptique, il existe ε > 0 petit tel que pour tout h ∈ E
vérifiant khk 6 ε, 21 d2 f (x∗ ).(h, h) + o(khk2 ) > α2 khk2 : autrement dit, pour khk assez petit on peut
absorber le terme en o. On obtient alors
α
f (x∗ + h) > f (x∗ ) + khk2
2
pour tout h ∈ E tel que khk 6 ε. Donc x∗ est un minimiseur local strict de f .
En fait, dans la condition suffisante, on a obtenu un résultat plus fort : non seulement le
minimiseur local est strict, mais de plus, la différence f (x∗ + h) − f (x∗ ) est au moins quadratique.
Remarque 16. Rappelons que, lorsque E = IRn , la Hessienne de f en x∗ est la matrice symétrique
2
réelle Hf (x∗ ) (carrée d’ordre n) formée par les dérivées partielles ∂x∂i ∂x
f
j
(x∗ ), i, j = 1, . . . , n. On
∗
rappelle que Hf (x ) est positive si
y > Hf (x∗ )y > 0 ∀y ∈ IRn ,
et Hf (x∗ ) est définie positive si elle est positive et de plus
y > Hf (x∗ )y = 0 ⇒ y = 0.
En dimension finie, il est toujours vrai que si Hf (x∗ ) est définie positive alors elle est α-elliptique,
pour un α > 0. En effet, la matrice Hf (x∗ ) est symétrique réelle, donc elle est diagonalisable en
base orthonormée. En notant λ1 > · · · > λn ses valeurs propres (qui sont réelles), on obtient donc
∀y ∈ IRn y > Hf (x∗ )y > λn kyk2 .
Ainsi, si Hf (x∗ ) est définie positive, alors λn > 0 et on conclut que Hf (x∗ ) est α-elliptique avec
α = λn = min Spec(Hf (x∗ )).
En dimension finie, la condition suffisante du deuxième ordre peut donc s’écrire : si ∇f (x∗ ) = 0
et si Hf (x∗ ) > 0 alors x∗ est un minimiseur local strict de f .
Mais l’implication “définie positive ⇒ α-elliptique" n’est pas vraie en dimension infinie, et c’est
pourquoi, dans le théorème, on a supposé que d2 f (x∗ ) est α-elliptique. En effet, il suffit d’imaginer
une Hessienne de taille infinie, et telle que, même si elle est positive, ses valeurs propres forment
une suite de réels strictement positifs décroissant vers 0. La conclusion du théorème peut alors être
mise en défaut.
Remarque 17. Même en dimension finie, si la Hessienne est seulement positive (mais pas α-
elliptique), la conclusion de minimiseur local est fausse. Par exemple, la fonction f : IR2 → IR
définie par
f (x, y) = x2 − y 4
est telle que
0 1 0
∇f (0, 0) = , Hf (0, 0) =
0 0 0
(donc, la Hessienne est positive mais pas définie positive) et (0, 0) est un point extremum qui
n’est ni un minimiseur ni un maximiseur. En effet, f (x, 0) = x2 donc dans la direction x on a un
minimum à l’origine. Mais par ailleurs, f (0, y) = −y 4 donc dans la direction y on a un maximum
à l’origine.
29
Remarque 18. Par le théorème 8, la condition d’α-ellipticité de la Hessienne en x∗ implique en
fait que, dans un voisinage de x∗ , la fonction f est α-convexe. C’est la raison pour laquelle on
trouve que non seulement x∗ est un minimiseur local strict, mais que, de plus, f (x∗ + h) − f (x∗ )
est au moins quadratique en h (pour khk petit).
Remarque 19. Les conditions obtenues ci-dessus sont locales, mais bien entendu, elles deviennent
globales si f est de plus convexe sur E.
Exemple 6. Un exemple important, déjà vu, est le cas de la fonction f : IRn → IR définie par
1 >
f (x) = x Ax − b> x
2
où A ∈ Mn (IR) est une matrice symétrique réelle définie positive et b ∈ IRn . L’unique minimiseur
x∗ ∈ IRn est déterminé par la condition ∇f (x∗ ) = 0, i.e., Ax∗ = b, donc, x∗ = A−1 b.
où A ∈ Mn,p (IR) est une matrice réelle ayant n lignes et p colonnes, et b ∈ IRn . Dans toute cette
section, k k est la norme euclidienne de IRn .
Il s’agit d’un problème de minimisation sans contrainte, pour la fonction f : IRp → IR définie
par
f (x) = kAx − bk2 .
La fonction f est de classe C ∞ , et en développant f (x + h) = f (x) + 2hAx − b, Ahi + kAhk2 , on
trouve
∇f (x) = 2A> (Ax − b) et Hf (x) = 2A> A.
On voit aussi que f est convexe puisque Hf (x) > 0 (en effet, h> Hf (x)h = 2kAhk2 ), mais elle n’est
pas forcément strictement convexe. La fonction f n’est pas forcément infinie à l’infini (en effet on
peut avoir ker A 6= {0}). Malgré cela, il est remarquable que le problème des moindres carrés a
toujours au moins une solution, comme l’affirme le résultat ci-dessous.
Théorème 18. Le problème des moindres carrés admet au moins une solution. Les solutions du
problème des moindres carrés sont exactement les solutions de “l’équation normale"
A> Ax = A> b.
De plus, si ker A = {0} (i.e., si A est injective) alors il existe une unique solution x ∈ IRp , qui est
30
Démonstration. Comme f est convexe, f a un minimiseur x si et seulement si ∇f (x) = 0, i.e.,
A> Ax = A> b. Cela montre que les solutions du problème des moindres carrés, si elles existent,
sont exactement les solutions de l’équation normale.
Notons que ker A = {0} ⇔ A> A inversible ⇔ A> A > 0 (en effet, il suffit de remarquer que
x A Ax = kAxk2 pour tout x ∈ IRp ). Sous cette condition, on a Hf (x) = 2A> A > 0, donc f est
> >
On a A# ∈ Mp,n (IR). Lorsque A est carrée inversible, on a A# = A−1 . Dans les autres cas, en
supposant A injective, la pseudo-inverse généralise l’inverse de A. Le cas typique est lorsqu’on
cherche résoudre un système de n équations à p inconnues
Ax = b
avec n > p (plus d’équations que d’inconnues) : en général, un tel système n’a pas de solution. En
pratique, il y a quantité de problèmes qui n’admettent pas de solutions (car on impose trop de
contraintes et on n’a pas assez de degrés de liberté), mais on cherche tout de même “le meilleur
compromis possible", i.e., on cherche à résoudre
Ax ' b
31
Lorsque ker A = {0}, on peut résoudre l’équation normale AA> x = A> b par la méthode de
Cholesky, la matrice AA> étant symétrique définie positive (et donc, elle s’écrit sous la forme LL>
avec L triangulaire inférieure).
Une autre méthode d’analyse numérique matricielle est la méthode QR : la décomposition
QR de la matrice A consiste à écrire A = QR où Q ∈ Mn (IR) est une matrice orthogonale et
R ∈ Mn,p (IR) est une matrice triangulaire supérieure. Le problème des moindres carrés est alors
équivalent à
minp kRx − Q> bk2
x∈IR
Mais, certainement, la méthode la plus efficace est de faire une décomposition SVD de la matrice
A. Etant donné l’importance de la décomposition SVD, et du fait qu’elle n’est hélas généralement
pas traitée en licence, cela vaut le coup d’y consacrer une section.
Remarque 23. Comme on l’a déjà remarqué, on a rg(A) = r, le nombre de valeurs singulières
non nulles. Par ailleurs, on a (d’après la remarque 21)
1/σ1 0
..
. 0
A # = U Σ# V > avec Σ# = 0 1/σr ∈ Mn,p (IR)
0 0
32
Remarque 24. Comme A = V ΣU > , premièrement, on a A> A = U Σ> ΣU > (avec U > = U −1 ),
donc la matrice U est exactement la matrice de passage qui diagonalise la matrice symétrique
réelle A> A, elle est donc composée par les vecteurs propres de A> A. La matrice carrée Σ> Σ est
diagonale, et les éléments de la diagonale sont exactement les valeurs propres de A> A (carrés des
valeurs singulières).
Deuxièmement, on a AA> = V ΣΣ> V > (avec V > = V −1 ), donc la matrice V est exactement la
matrice de passage qui diagonalise la matrice symétrique réelle AA> , elle est donc composée par
les vecteurs propres de AA> . La matrice carrée ΣΣ> est diagonale, et les éléments de la diagonale
sont exactement les valeurs propres de AA> (carrés des valeurs singulières) : rien d’étonnant car
en fait les valeurs propres non nulles de A> A et de AA> sont les mêmes.
Cette remarque donne presque la preuve du théorème.
Démonstration. Montrons le théorème pour n > p (pour n 6 p, il suffit d’appliquer la décomposi-
tion SVD à A> ). On note ui les vecteurs propres de la matrice symétrique réelle A> A, associés aux
valeurs propres σi2 , i = 1, . . . , p. On définit la matrice orthogonale U ∈ Mp (IR) dont les colonnes
sont les ui . On a alors exactement A> A = U Σ> ΣU > .
En particulier, cette égalité donne hAuj , Aui i = u> > 2
j A Aui = σi δi,j (avec δi,j = 1 si i = j
et 0 sinon) pour i = 1, . . . , p ; notons que Aui = 0 pour i = r + 1, . . . , p (car σi = 0). Donc,
en posant vi = σ1i Aui pour i = 1, . . . , r, la famille (v1 , . . . , vr ) est orthonormée. On la complète
en une base orthonormée (v1 , . . . , vn ) de IRn (de manière quelconque) et on définit la matrice
orthogonale V ∈ Mn (IR) dont les colonnes sont les vi . On a Aui = σi vi pour i = 1, . . . , p. Cela
donne exactement AU = V Σ.
33
k 6 min(n, p), cherchons une matrice Ak de rang k qui soit une “bonne approximation" de la
matrice A. La matrice
σ1 0
..
. 0
A k = V Σk U > avec Σk = 0 σk ∈ Mn,p (IR)
0 0
est la meilleure approximation de A parmi les matrices de rang k (on n’a pas forcément unicité)
au sens de la norme k k2 (qui est la norme matricielle subordonnée à la norme euclidienne), et
on a kA − Ak k2 = σk+1 . Si A est la matrice d’une image, son approximation de rang k est une
compression de l’image. Amusez-vous à coder cela (en Matlab, Scilab ou Python), et prenez diverses
valeurs de k. Vous constaterez que, même pour k petit on arrive à avoir une image de bonne qualité.
Avant l’invention des normes jpeg, la décomposition SVD a longtemps servi pour la compression
d’image. Elle est désormais supplantée, en efficacité, par les ondelettes.
P (ti ) ' bi , i = 1, . . . , n.
Bien entendu, lorsque n 6 p, on peut trouver un polynôme qui réalise exactement toutes les égalités
(on utilise la théorie de l’interpolation de Lagrange). Mais, ici, on a plutôt p n et les égalités
ne peuvent pas toutes être satisfaites. On cherche donc un meilleur compromis, en cherchant le
polynôme P ∈ IRp−1 [X] qui minimise
n
X
f (P ) = (P (ti ) − bi )2 .
i=1
Il s’agit d’un problème de moindres carrés. En effet, soit (φj )06j6p−1 une base de IRp−1 [X] (on
n’est pas obligé de considérer la base canonique). On cherche
p−1
X
P = aj φj .
j=0
En posant
a0 b1
x = ... ∈ IRp , b = ... ∈ IRn ,
A = φj (ti ) ∈ Mn,p (IR)
16i6n
06j6p−1
ap−1 bn
34
La matrice A est injective car les ti sont tous distincts. Ce problème de moindres carrés admet
donc une unique solution, qui est x = A# b.
Exemple 7. Prenons p = 2, φ0 = 1 et φ1 = X (base canonique), et retrouvons les formules
connues de régression linéaire. On a ici
Xn X n
1 t1 n ti b i
1 ··· 1
.. .. , A> = > i=1
A> b = X
i=1
A = . , A A = X ,
. t1 ··· tn n Xn n
2
1 tn ti ti ti bi
i=1 i=1 i=1
a0
et comme x = = (AA> )−1 A> b, on trouve
a1
n
X n
X n
X n
X n
X n
X n
X
t2i bi − ti ti bi − ti bi + n t i bi
i=1 i=1 i=1 i=1 i=1 i=1 i=1
a0 = !2 , a1 = !2
n
X n
X n
X n
X
n t2i − ti n t2i − ti
i=1 i=1 i=1 i=1
Cela nous amène naturellement à choisir ρk > 0 petit et une direction de descente dk ∈ IRn vérifiant
35
par exemple, dk = −∇f (xk ). Avec un tel choix, si ρk est assez petit alors on a bien f (xk+1 ) < f (xk ),
et on peut alors espérer que, lorsque k → +∞, on ait xk → x∗ , à condition, bien sûr, d’avoir choisi
x0 pas trop loin de x∗ (car, si f n’est pas convexe, elle peut avoir d’autres minimiseurs locaux).
Une façon de choisir le pas ρk > 0 est de minimiser la fonction ϕ(t) = f (xk + tdk ) :
ρk = argmin ϕ(t)
t>0
Un tel minimum existe par exemple si f est infinie à l’infini : en effet, d’une part, ϕ0 (0) =
df (xk ).dk = h∇f (xk ), dk i < 0, donc ϕ décroît à partir de t = 0, au moins sur un petit inter-
valle [0, ε], et d’autre part, ϕ(t) → +∞ lorsque t → +∞, donc ϕ a un minimum. En un tel
minimiseur ρk , on doit avoir
ce qui nous conduira, un peu plus loin, à définir la méthode de gradient à pas optimal.
Sans chercher, pour le moment, un pas optimal, on peut au moins choisir un pas ρk = ρ > 0
fixe, assez petit. Cela donne la méthode du gradient à pas fixe, qu’on va analyser ci-après. On peut
aussi choisir des pas variables.
On peut choisir des directions de descente autres que dk = −∇f (xk ). Par exemple, dans la
méthode de gradient conjugué, on construira une suite de directions de descente, ayant certaines
propriétés d’orthogonalité. Cette construction astucieuse assure une convergence en nombre fini
d’itérations pour des fonctions quadratiques.
Les variantes sont nombreuses, on va en voir et en analyser quelques-unes.
De manière générale, l’algorithme d’une méthode de descente, ainsi que définie ci-dessus, est le
suivant :
1. Initialisation, k = 0 : choisir x0 ∈ IRn et ε > 0.
2. Itération k : ayant choisi un pas ρk et une direction de descente dk , on calcule
xk+1 = xk + ρk dk .
k∇f (xk )k 6 ε
kxk+1 − xk k 6 ε
c’est-à-dire, lorsque les itérés ne progressent plus suffisamment ; ou bien (variante), lorsque
kf (xk+1 ) − f (xk )k 6 ε
36
Théorème 20. Soit U un ouvert convexe borné de IRn sur lequel f a un unique minimiseur x∗ . On
suppose que, sur U , la fonction f est de classe C 1 , strictement convexe 1 , et que f est à gradient
Lipschitz 2 , i.e., qu’il existe M > 0 tel que
∀x, y ∈ U k∇f (x) − ∇f (y)k 6 M kx − yk. (3.1)
Soient 0 < β1 < β2 < 2
[Link] η > 0 tel que V = {x ∈ U | f (x) 6 f (x∗ ) + η} ⊂ U . Si, pour
tout k ∈ IN, le pas ρk vérifie ρk ∈ [β1 , β2 ], alors la méthode de gradient, initialisée à un x0 ∈ V
arbitraire, converge vers x∗ : la suite (xk )k∈IN définie par
xk+1 = xk − ρk ∇f (xk ), x0 ∈ V
reste dans V et converge vers x∗ .
On appelle cette méthode la méthode de gradient à pas variable. Le choix du pas à chaque
itération peut être fait selon différents critères (on en verra plus loin). Lorsque le pas ρk = ρ est
fixe, on parle de méthode de gradient à pas fixe.
Bien entendu, en pratique on ne connait pas à l’avance un ouvert convexe U sur lequel f
aurait un unique minimiseur. La recherche globale de telles “régions favorables" est souvent faite
de manière heuristique (à moins évidemment que f est globalement convexe). Souvent, également,
on reste un gradient à pas fixe, on choisit une initialisation un peu au hasard et on espère que
l’algorithme va converger. Mais, lorsqu’on a plus d’informations sur la fonction, il faut essayer de
s’en servir !
Remarque 26. On a vu à la remarque 7 (dont on remet les figures ici) que, comme f est supposée
convexe sur U , les ensembles de sous-niveau {x ∈ U | f (x) 6 α} sont convexes. De plus, si f est
C 1 en x et si f (x) = α, le gradient ∇f (x) est un vecteur qui est orthogonal à l’ensemble de niveau
f = α (qui est une hypersurface au voisinage de x). Cette interprétation du gradient est à la base
de la méthode de gradient : le gradient est orthogonal aux surfaces de niveau, et est orienté dans
le sens où f croît, comme on le voit sur la figure.
Sur la figure de droite, les ensembles de niveau f = α sont en bleu (les ensembles de sous-
niveau sont leur intérieur), au centre se trouve le minimiseur x∗ . En rouge sont représentées les
courbes intégrales du gradient de f . Ce sont les courbes qui sont, en approximation, suivies par les
itérations xk .
L’ensemble V du théorème est un ensemble de sous-niveau compact (fermé borné) de f , contenu
dans U (et contenant le minimiseur x∗ ), il est donc convexe fermé. Son intérieur V̊ = {x ∈ U |
f (x) < f (x∗ ) + η} est un voisinage ouvert convexe de x∗ .
1. Cela est vrai si f est C 2 et Hf (x) > αIn pour tout x ∈ U , pour un α > 0.
2. Cela est vrai si Hf (x) 6 M In pour tout x ∈ U .
37
Démonstration. L’unique minimiseur x∗ de f sur U est caractérisé par ∇f (x∗ ) = 0 (par stricte
convexité). L’ensemble V est convexe compact et contient x∗ en son intérieur. Supposons être
à l’itération k, avec xk ∈ V (c’est le cas pour k = 0, ce qui initialise la récurrence). On pose
xk+1 = xk − ρk ∇f (xk ), et on suppose que ∇f (xk ) 6= 0 (sinon il n’y a rien à faire). Attention, on
ne sait pas encore que xk+1 ∈ V : en effet, il se pourrait que ρk > 0 soit trop grand et que xk+1
sorte du convexe V . On va voir que ce n’est pas le cas grâce au choix de ρk , mais pour prendre en
compte cette difficulté, on pose
et bien sûr on a x(ρk ) = xk+1 . Notons que, comme ∇f (xk ) 6= 0, on a x(ρ) 6= xk pour tout ρ > 0. Il
existe ε > 0 assez petit, tel que, pour tout ρ ∈]0, ε], on ait x(ρ) ∈ V (voir figure dans la remarque
ci-dessus) : pour le montrer rigoureusement, il suffit d’écrire un développement limité au premier
ordre, f (x(ρ)) = f (xk − ρ∇f (xk )) = f (xk ) − ρk∇f (xk )k2 + o(ρ) < f (xk ) pour ρ > 0 assez petit,
et donc x(ρ) ∈ V , et même, x(ρ) ∈ V̊ . On pose
(le supremum est bien un maximum car V est fermé). On vient de voir que β > 0, et notre objectif
est de montrer que β = β2 .
Soit ρ ∈]0, β] arbitraire. Suivant la technique déjà largement utilisée, posons
ce qui correspond à regarder la fonction f le long du segment [xk , x(ρ)] : on a ϕ(0) = f (xk ),
R1
ϕ(1) = f (x(ρ)) et ϕ0 (t) = h∇f (xk + t(x(ρ) − xk )), x(ρ) − xk i. On a ϕ(1) = ϕ(0) + 0 ϕ0 (t) dt, ce
qui donne
Z 1
f (x(ρ)) = f (xk ) + h∇f (xk + t(x(ρ) − xk )), x(ρ) − xk i dt
0
Z 1
= f (xk ) + h∇f (xk ), x(ρ) − xk i + h∇f (xk + t(x(ρ) − xk )) − ∇f (xk ), x(ρ) − xk i dt
0
et comme ρ 6 β2 , on a
M 1
f (x(ρ)) − f (xk ) 6 − kx(ρ) − xk k2 < 0
2 β2 | {z }
| {z } ρ2 k∇f (xk )k2 >0
<0
38
Ainsi, tout le raisonnement ci-dessus peut être appliqué à ρ = ρk , puisque, par hypothèse,
ρk 6 β2 . On a donc
M 1
f (xk+1 ) − f (xk ) 6 − kxk+1 − xk k2 < 0
2 β2 | {z }
| {z } ρ2 k∇f (xk )k2 >0
k
<0
dont on déduit deux choses. Premièrement, la suite (f (xk ))k∈IN est décroissante. Comme elle est
minorée (f a un minimum sur U ), elle converge. Donc la suite de réels positifs f (xk ) − f (xk+1 )
converge vers 0. Deuxièmement, en renversant l’inégalité, on a
1
kxk+1 − xk k2 6 1
(f (xk ) − f (xk+1 ))
β2−M 2
| {z }
>0
Remarque 27. Si dans le théorème on suppose que f est strictement convexe sur IRn tout entier,
que f est infinie à l’infini, et que l’hypothèse (3.1) est globale sur IRn , alors on n’a pas besoin de
faire le raisonnement avec x(ρ) : on peut l’appliquer directement à xk+1 .
Le théorème 20 est plus fort et s’applique aux minimiseurs locaux. Il montre aussi que, pour
assurer la convergence de la méthode de gradient, on a intérêt à avoir une idée a priori d’où se
situent les minimiseurs locaux (et en effet, en pratique, on a souvent une connaissance intuitive du
problème qui permet de faire une première localisation grossière) et à initialiser la méthode en un
point x0 en lequel on a la propriété locale de convexité.
Remarque 28 (Méthode du gradient à pas fixe). Lorsque ρk = ρ est constant, on parle de méthode
de gradient à pas fixe. C’est certainement la méthode d’optimisation sans contrainte la plus simple.
Il est intéressant de noter que la méthode de gradient à pas fixe s’interprète comme une méthode
de point fixe. En effet, en posant
Fρ (x) = x − ρ∇f (x),
on constate que
xk+1 = xk − ρ∇f (xk ) ⇔ xk+1 = Fρ (xk )
Autrement dit, les itérations de la méthode de gradient à pas fixe sont exactement les itérations
de la méthode du point fixe de Picard.
Or, on sait que la méthode du point fixe converge lorsque l’application Fρ est K-contractante,
i.e., K-Lipschitzienne avec 0 6 K < 1. Vérifions que c’est bien le cas pour Fρ , sous la condition
(3.1) du théorème 20 (i.e., f est à gradient Lipschitz) et sous la condition, plus forte que dans le
39
théorème, que f soit α-convexe pour un α > 0 : on a, pour tous x, y ∈ V ,
6 (1 − 2ρα + ρ2 M 2 ) kx − yk2
| {z }
K2
2
2α α
On a K < 1 si et seulement si 0 < ρ < M 2 . De plus, la valeur minimale de K est K = 1 − M 2 et
α
est obtenue pour ρ = M 2 .
Or, on sait que, dans la méthode de point fixe, on a xk+1 − x∗ = Fρ (xk ) − Fρ (x∗ ) (car x∗ =
Fρ (x∗ )), donc kxk+1 − x∗ k 6 Kkxk − x∗ k et donc
kxk − x∗ k 6 K k kx0 − x∗ k ∀k ∈ IN
ce qui veut dire que la méthode de gradient à pas fixe converge linéairement, avec vitesse de
convergence K (voir le rappel ci-dessous).
Ainsi, dans le cas d’un gradient à pas fixe, si f est α-convexe et à gradient Lipschitz sur U ,
on a obtenu (certes, sous des hypothèses un peu plus fortes que dans le théorème 20) une preuve
alternative, qui a l’avantage de donner unep vitesse de convergence explicite : la méthode converge
2α
linéairement si ρ < M 2 , à vitesse K = 1 − 2ρα + ρ2 M 2 . De plus, le meilleur pas fixe possible
(celui qui donne le K le plus petit) est ρ = Mα2 .
2α
Notons que α < M , et donc le seuil de convergence M 2 établi ici est bel et bien inférieur au
2
seuil de convergence M établi dans le théorème 20.
Remarque 29 (Vitesse de convergence). On rappelle le vocabulaire suivant. Soit (xk )k∈IN une
suite convergeant vers x∗ (on suppose ci-dessous que xk 6= x∗ pour tout k).
On dit que (xk )k∈IN converge à l’ordre p > 1 vers x∗ , à vitesse c ∈]0, 1[, si
kxk+1 − x∗ k
lim = c. (3.2)
k→+∞ kxk − x∗ kp
xk+1 = xk + ρk dk
où ρk = argmin f (xk + ρdk )
dk = −∇f (xk ) ρ>0
40
Posons
ϕ(ρ) = f (xk + ρdk ).
Il s’agit donc, à chaque étape, de déterminer (s’il existe) le minimum de la fonction ϕ sur ]0, +∞[.
Notons qu’un tel minimum existe si f est infinie à l’infini, car ϕ0 (0) = −kdk k2 < 0 (pourvu que
dk = ∇f (xk ) 6= 0 bien sûr) donc ϕ décroît pour ρ > 0 petit, et par ailleurs ϕ(ρ) → +∞ lorsque
ρ → +∞. On a ϕ0 (ρk ) = 0.
Déterminer un tel minimum global peut toutefois s’avérer coûteux, et en pratique, on cherche
donc une approximation du pas optimal ρk > 0. En supposant que ρk est petit, on écrit le déve-
loppement limité
kdk k2
ρk = >
où Hk = Hf (xk ).
dk Hk dk
Notons que ces formules sont exactes lorsque f est quadratique, i.e., lorsque f (x) = 12 x> Ax − b> x
avec A symétrique définie positive.
Même dans le cas non linéaire, on a l’habitude d’appeler méthode de gradient à pas optimal la
méthode suivante :
Cette méthode (dont on ne fait pas l’analyse de convergence ici) est excellente en pratique mais
peut toutefois s’avérer coûteuse en grande dimension car elle réclame de calculer la Hessienne de
f à chaque étape. Or, dans des problèmes de sciences des données, cela peut être rhédibitoire. Il
faut donc chercher d’autres alternatives.
On verra, dans les méthodes de type Newton, qu’on pourrait remplacer la Hessienne par une
“quasi-Hessienne" (i.e., on remplace Hk par une suite dont on sait démontrer qu’elle constitue,
pour k grand, une bonne approximation de la Hessienne).
Mais, même le temps de calcul du gradient peut être problématique lorsque le nombre de
variables est très grand. Certains problèmes comportent des millions de variables, et si la fonction
est un peu trop non linéaire, calculer le gradient complet à chaque étape est trop lent : en effet, si
n = 106 par exemple, chaque itération k réclame de calculer n = 106 dérivées partielles (évaluées
au point xk ). Il faut alors trouver des moyens de réduire la dimension, d’une manière ou d’une
autre. Dans la section suivante, on décrit une méthode très simple de gradient coordonnée par
coordonnée, ou bloc par bloc, qui permet, en fait, “d’avancer de manière un peu diagonale".
qu’on appelle méthode de relaxations successives de type Jacobi, ou, comme variante plus utilisée :
min f (xk+1
1 , . . . , xk+1 k k
i−1 , x, xi+1 , . . . , xn ), pour i = 1, . . . , n
x∈IR
41
qu’on appelle méthode de relaxations successives de type Gauss-Seidel. Ici, par commodité, l’indice
d’itération k est mis en haut ; l’indice du bas étant la coordonnée.
Ces dénominations Jacobi / Gauss-Seidel proviennent des méthodes itératives de résolution
de systèmes linéaires Ax = b, qui consistent à décomposer A = D − E − F avec D diagonale, E
triangulaire inférieure et F triangulaire supérieure (qu’on appelle de manière générale les méthodes
de relaxation). La méthode de Gauss-Seidel est généralement préférée – plus efficace car, à l’étape
k, au fur et à mesure des calculs sur i = 1, . . . , n, on se sert des valeurs de xk+1
i déjà calculées.
Quoi qu’il en soit, l’idée est ici, à chaque étape k, de remplacer l’itération “globale sur x" (qui
consiste, dans une méthode de descente, à calculer xk+1 = xk + ρk dk où dk est une direction de
descente, par exemple dk = −∇f (xk )) par n calculs qui se font chacun, coordonnée par coordonnée.
Par exemple, à l’étape k, pour chaque i = 1, . . . , n, on peut choisir une méthode de gradient (à pas
ρk ), ce qui donne, dans la méthode de Gauss-Seidel :
∂f k+1
xk+1 = xki − ρk (x , . . . , xk+1 k k k
i−1 , xi , xi+1 , . . . , xn ), i = 1, . . . , n
i
∂xi 1
Notons que, par rapport à une méthode de gradient classique, on calcule le même nombre de
dérivées partielles ! Mais les calculs et itérations ne sont pas faits dans le même ordre, et cela
change tout. Dans la méthode de gradient classique, à chaque étape k, on calcule toutes les dérivées
partielles pour i = 1, . . . , n et on fait un pas de gradient. Alors que, dans la méthode coordonnée
par coordonnée, à chaque étape k, on itère sur i = 1, . . . , n, et pour chaque i, on fait un pas de
gradient par rapport à la coordonnée xi . C’est différent ! On conçoit en effet que cela peut aller
plus vite car, à chaque étape, on peut espérer avoir significativement avancé dans la convergence
par rapport à la coordonnée xi .
Au lieu de choisir une méthode de gradient, on peut bien entendu choisir toute autre méthode
de minimisation de fonctions à une variable, ce qui engendre plein de variantes possibles.
D’autre part, ci-dessus, on s’est ramené à des minimisations de fonctions à une variable, mais on
aurait pu aussi faire des regroupements de variables, des “blocs", et remplacer l’itération i = 1, . . . , n
qui est faite sur toutes les coordonnées, par une itération faite sur les blocs de variables. On appelle
ce type de méthode, les méthodes de descente bloc par bloc. En fait la méthode s’écrit exactement
de la même façon que ci-dessus pourvu de remplacer xi ∈ IR par xi ∈ IRni où ni est la dimension
du bloc i.
Théorème 21. On suppose que f : IRn → IR est C 1 , infinie à l’infini et strictement convexe. Alors
la méthode de gradient coordonnée par coordonnée de type Gauss-Seidel converge vers l’unique
minimiseur de f .
La preuve (qu’on peut trouver, par exemple, dans le livre de Glowinski, Lions, Trémolières,
1981, Chapter 2) est du même type que celle du théorème 20.
Enfin, dernière remarque, ces méthodes ont connu un fort regain d’intérêt récemment avec
l’avènement du deep learning. Le choix des coordonnées ou des blocs de descente peut être fait de
manière aléatoire, par exemple avec probabilité uniforme. La convergence de ce type de méthode
a été étudiée récemment.
42
Dans cette section, on discute de méthodes heuristiques permettant d’orienter le choix du pas
ρk à chaque étape. On a déjà vu qu’une (très) bonne méthode pour choisir le pas ρk était la
méthode du gradient à pas optimal. Mais elle présente le défaut de devoir calculer la Hessienne
(ou au moins une approximation de la Hessienne) à chaque itération, ce qui peut être prohibitif si
n est très grand.
On va donc décrire ici des possibilités très simples qui permettent d’orienter le choix du pas ρk ,
dans une méthode de descente xk+1 = xk + ρk dk , de manière intelligente. Comme précédemment,
à l’étape k, on pose
ϕ(ρ) = f (xk + ρdk ).
L’objectif d’une recherche linéaire est de ne pas chercher à calculer le minimum de ϕ (car cela peut
être trop coûteux), et de déterminer un pas ρk permettant d’assurer que ϕ décroît suffisamment.
En même temps, il faut assurer que ρk ne soit ni trop petit ni trop grand. Tout cela est très
heuristique ! Et largement basé sur l’intuition.
Dans les différentes règles décrites ci-dessous, on décide de choisir ρk dans un intervalle [a, b]
appelé “intervalle de sécurité" : autrement dit, un réel ρ < a est considéré comme étant trop petit,
un réel ρ > b est considéré comme étant trop grand, et lorsqu’on trouve un ρ ∈ [a, b] on décide
qu’il est convenable.
Règle d’Armijo. On fixe un réel m ∈]0, 1[, et on décide que ρ convient si ϕ(ρ) 6 ϕ(0)+mϕ0 (0)ρ.
Autrement dit, ρ convient si ϕ(ρ) est en dessous de la droite de pente m ϕ0 (0) passant par
(0, ϕ(0) = f (xk )) (rappelons que ϕ0 (0) < 0 grâce au choix de la direction de descente).
Dans la règle d’Armijo, on a donc a = 0, et b est l’abscisse du premier point d’intersection entre
la droite de pente m ϕ0 (0) ci-dessus et le graphe de ϕ. Mais comme a = 0, on n’a pas de borne
inférieure sur le pas, et donc, souvent, on combine cette règle à d’autres règles.
Règle de Goldstein. On fixe deux réels 0 < m1 < m2 < 1, et on décide que ρ convient si
autrement dit, si ϕ(ρ) est en dessous de la droite de pente m1 ϕ0 (0) passant par (0, ϕ(0) = f (xk )),
et au-dessus de la droite de pente m2 ϕ0 (0) passant par (0, ϕ(0) = f (xk )).
43
On choisit en général m2 > 12 de façon à ce que le pas optimal déterminé dans le cas quadratique
appartienne à l’intervalle de sécurité (vérifiez cela par un calcul). Par exemple, m1 = 0.1 et m2 =
0.7.
Règle de Wolfe. On fixe deux réels 0 < m1 < m2 < 1, et on décide que ρ convient si
autrement dit, si ϕ(ρ) est en dessous de la droite de pente m1 ϕ0 (0) passant par (0, ϕ(0) = f (xk )),
et si la pente de ϕ en 0 est plus grande que m2 ϕ0 (0).
Dans cette règle, il faut estimer ϕ0 (0) et donc calculer ∇f (xk ) à chaque pas. Mais, si cela n’est
pas trop coûteux, cette information différentielle est plus précise.
Bien sûr, on peut imaginer quantité d’autres variantes.
F (x) = ∇f (x)
F (x) = 0 (3.3)
44
qui est un système (non linéaire) de n équations à n inconnues.
Bien entendu, ce système étant non linéaire, il peut admettre plusieurs solutions. Nous sommes
intéressés à trouver la solution x∗ , où x∗ est un minimiseur (au moins local) de f . Mais notons que
tout minimiseur local, tout maximiseur local vérifie (3.3) – et pas seulement : on peut aussi avoir
des solutions de (3.3) qui ne sont ni un minimiseur ni un maximiseur local (comme par exemple 0
pour l’équation x3 = 0).
Le problème (3.3) ne peut donc être “bien posé" (i.e., admettre une unique solution) que dans
un ouvert U de IRn , éventuellement assez petit. Et encore, malgré ce caractère bien posé local, le
fait d’être l’unique solution dans U de (3.3) ne distinguera pas entre un minimum, maximum ou
extremum qui n’est aucun des deux.
ce qui nous conduit à choisir xk+1 tel que F (xk ) + dF (xk ).(xk+1 − xk ) = 0, i.e., lorsque c’est
possible,
xk+1 = xk − dF (xk )−1 .F (xk )
Ici, dF (xk ), la différentielle de F au point xk , s’identifie à une matrice carrée de taille n, et on doit
donc assurer qu’elle soit inversible. Comme xk est supposé proche de x∗ , cela conduit naturellement
à supposer que dF (x∗ ) est inversible. Cette condition est l’une des conditions suffisantes assurant
le caractère localement bien posé de la méthode de Newton.
Théorème 22. On suppose que F : IRn → IRn est de classe C 1 , que F (x∗ ) = 0, que dF (x∗ ) est
inversible et que dF est Lipschitzienne 3 au voisinage de x∗ , i.e., il existe r0 > 0 et L > 0 tel que
Alors il existe r ∈]0, r0 ] tel que, pour tout x0 ∈ B(x∗ , r), la suite (xk )k∈IN définie par l’itération
partant de x0 ∈ B(x∗ , r), est bien définie, reste dans la boule B(x∗ , r), et converge vers x∗ . De
plus, la convergence est quadratique (voir remarque 29), i.e., il existe c > 0 tel que
kxk+1 − x∗ k 6 ckxk − x∗ k2 ∀k ∈ IN
Dans (3.4), k k est (par exemple) la norme Euclidienne de IRn , la boule B(x∗ , r0 ) est la boule
Euclidienne de centre x∗ et de rayon r0 . La norme k kL(IRn ) est la norme d’opérateur, dont on
rappelle qu’elle est définie par
k`(x)k
k`kL(IRn ) = sup = sup k`(x)k
x∈IRn \{0} kxk x∈IRn
kxk=1
45
Démonstration. Comme F est C 1 et dF (x∗ ) est inversible, il existe 0 < R < r0 tel que dF (x) reste
inversible pour tout x ∈ B(x∗ , R). On pose alors
Soit x0 ∈ B(x∗ , R). Par définition, on a x1 = x0 − dF (x0 )−1 .F (x0 ), et comme F (x∗ ) = 0 on peut
écrire
x1 − x∗ = x0 − x∗ − dF (x0 )−1 .(F (x0 ) − F (x∗ ))
= −dF (x0 )−1 .(F (x0 ) − F (x∗ ) − dF (x0 ).(x0 − x∗ ))
R1
Mais d’une part on a F (x0 ) − F (x∗ ) = 0 dF (x∗ + t(x0 − x∗ )).(x0 − x∗ ) dt (formule de Taylor avec
R1
reste intégral à l’ordre 1), et d’autre part, trivialement, dF (x0 ).(x0 − x∗ ) = 0 dF (x0 ).(x0 − x∗ ) dt,
donc
Z 1
F (x0 ) − F (x∗ ) − dF (x0 ).(x0 − x∗ ) = dF (x∗ + t(x0 − x∗ )) − dF (x0 ) .(x0 − x∗ ) dt
0
46
Remarque 33 (Rapport avec le point fixe de Banach). Résoudre F (x) = 0 est équivalent à
résoudre
GM (x) = x − M F (x) = x
où M ∈ GLn (IR) est une matrice inversible quelconque, autrement dit, à trouver un point fixe de
GM . L’algorithme du point fixe consiste à considérer la suite définie par récurrence
On sait que cette suite converge vers x∗ lorsque l’application GM est contractante. Or, dGM (x) =
In − M dF (x), et donc, pour assurer que GM soit contractante dans un voisinage de x∗ , on peut
choisir M = dF (x∗ )−1 ; mais comme on ne connaît pas x∗ , il est naturel de choisir M = Mk
dépendant de k, en posant Mk = dF (xk )−1 , ce qui est exactement la méthode de Newton !
Remarque 34. A moins que la fonction F vérifie certaines conditions globales, la convergence de
la méthode de Newton est locale. On a un domaine de convergence, qui est souvent assez petit.
La difficulté majeure de la méthode de Newton est donc son initialisation : il faut en effet être
capable de deviner un point initial x0 qui est dans le domaine de convergence (donc, proche de
x∗ ). Cela peut sembler paradoxal : pour déterminer x∗ qui est inconnu, on a intérêt à déjà en
connaître une approximation, si on veut assurer la convergence de la méthode de Newton. Ce
genre de considération est inhérent à toute méthode locale.
De nombreuses méthodes existent pour rendre la méthode de Newton “un peu moins locale",
et arriver à la faire converger même en n’ayant qu’une idée assez vague d’où se situe le point x∗
recherché. Par exemple :
— Il existe des méthodes de Newton “globales", qui requièrent toutefois des hypothèses fortes
sur la fonction F . On n’en parle pas ici.
— Avant d’appliquer une méthode de Newton, on peut, au préalable, appliquer une méthode de
descente à la fonction
x 7→ kF (x)k2 .
L’avantage est que cette méthode peut converger “plus facilement", avec une initialisation
grossière. Ainsi, en quelques itérations, elle peut fournir un point un peu meilleur dont on
peut se servir pour initialiser (on l’espère, avec succès) la méthode de Newton.
On appelle cela une méthode hybride. Il y a évidemment plein de manières d’hybridifier la
méthode de Newton. A vous de combiner astucieusement les méthodes que vous connaissez,
en fonction du problème !
47
— Une méthode très puissante pour faire converger la méthode de Newton est de la combiner à
une méthode de continuation : on déforme l’équation F (x) = 0 à l’aide d’un paramètre (qui,
souvent, est un ou plusieurs paramètres qui dans le problème rendent délicate sa résolution
numérique, et qu’il s’agit alors, d’une manière ou d’une autre, de relaxer). Appelons λ ce
paramètre, et supposons pour simplifier que λ ∈ [0, 1] (mais il pourrait y avoir plusieurs
paramètres !). On suppose maintenant que la fonction F dépend aussi du paramètre λ, i.e.,
on a F (λ, x), avec, pour λ = 1, F (1, x) = F (x) qui est la fonction de départ. On veut
maintenant résoudre
F (λ, x) = 0
dont l’unique solution locale (sous des conditions qui assurent le caractère bien posé) est
notée xλ . Pour λ = 1, on a x1 = x∗ qui est le minimiseur recherché. Pour λ = 0, on suppose
que l’équation
F (0, x) = 0
est “facile à résoudre" (c’est ce qu’on doit être capable d’assurer en choisissant adéquatement
le paramètre de continuation), sa solution étant x0 . C’est notre point de départ. A partir de
ce point de départ, on résout maintenant F (λ, x) = 0 pour λ petit (par exemple, λ = 0.1),
par une méthode de Newton initialisée à x0 qui était la solution de F (0, x0 ) = 0. Comme
λ est petit, on peut s’attendre à ce que xλ ' x0 , donc, à ce que la méthode converge. Si la
méthode a convergé, on recommence avec λ un peu plus grand (par exemple, λ = 0.2), et
ainsi de suite. En cas d’échec, on diminue l’incrément de λ.
On voit ainsi qu’au lieu de résoudre un seul problème de Newton, on en résout toute une série.
Comme l’exécution de la méthode de Newton est quasi-instantanée (vitesse de convergence
quadratique), cela n’est pas un problème et en général les continuations sur des problèmes
de Newton sont rapides et efficaces. Il en existe de nombreuses variantes.
Remarque 35. Appliqué à F (x) = ∇f (x), l’algorithme de Newton s’écrit
On voit qu’il est nécessaire, à chaque itération, de calculer la Hessienne de f en xk . Cela peut être
coûteux, surtout si on est en grande dimension. Cela nous conduit aux méthodes de quasi-Newton.
Ainsi, on souhaite déterminer une suite de matrices (H̃k )k∈IN ou (Bk )k∈IN telles que
yk = H̃k sk ou Bk yk = sk .
Mais il faut déterminer de telles matrices, qui soient symétriques définies positives, pour tout k !
48
Algorithme DFP (Davidon, Fletcher, Powell). Cette formule de mise à jour est une formule
de correction de rang 2 donnée par
sk s>
k Bk yk yk> Bk
Bk+1 = Bk + −
s>
k yk yk> Bk yk
On peut montrer (c’est admis, ici) que, pour toute matrice B0 symétrique définie positive (par
exemple, B0 = I), l’algorithme de quasi-Newton (dit DFP)
sk s>
k Bk yk y > Bk
xk+1 = xk − Bk ∇f (xk ), Bk+1 = Bk + >
− > k
sk yk yk Bk yk
(qui est l’algorithme de Newton dans lequel on a remplacé Hk−1 par Bk ), converge vers un minimum
local x∗ de f , sous les mêmes hypothèses que la méthode de Newton, et de plus,
Algorithme BFGS (Broyden, Fletcher, Goldfarb, Shanno). C’est aussi une formule de
correction de rang 2, qui consiste à intervertir les rôles de sk et yk dans la formule DFP. Elle fournit
alors une approximation H̃k de la Hessienne Hk . L’itération est définie par
On peut montrer (c’est admis, ici) que, pour toute matrice H̃0 symétrique définie positive (par
exemple, H̃0 = I), l’algorithme de quasi-Newton (dit BFGS)
(qui est l’algorithme de Newton dans lequel on a remplacé Hk par H̃k ), converge vers un minimum
local x∗ de f , sous les mêmes hypothèses que la méthode de Newton, et de plus,
Par rapport à la méthode DFP, la méthode BFGS nécessite de calculer l’inverse de la matrice H̃k ,
et peut donc sembler moins intéressante. Toutefois, elle s’avère être en général meilleure et plus
robuste que la méthode DFP. La méthode BFGS est la méthode de quasi-Newton la plus connue
et la plus utilisée.
49
Comme la fonction ρ 7→ f (xk − ρgk ) a, par définition, un minimum en ρk , sa dérivée doit être nulle
en ρk , donc ∇f (xk −ρk gk )> gk = 0, mais, par un développement limité à l’ordre 1, ∇f (xk −ρk gk ) =
∇f (xk ) − ρk Hf (xk )gk + o(ρk ), donc, en approximation, (gk − ρk Hk gk )> gk = 0 ce qui conduit à
prendre
kgk k2
ρk = >
gk Hk gk
qui est ce qu’on a appelé le “pas optimal".
Dans la méthode de Barzilai Borwein, on considère l’itération xk+1 = xk − ρk gk = xk − (ρk I)gk
où I est la matrice identité, et on voudrait (c’est sérieusement gonflé !) pouvoir choisir le pas ρk
tel que
ρk I ' Hk−1
car cela nous rapprocherait de la méthode de Newton, connue pour converger très vite. Bien sûr,
cette pseudo-égalité n’a aucun sens car la matrice Hk−1 n’est pas diagonale. Pourtant, un peu à
la manière des moindres carrés, on va chercher à faire en sorte qu’elle soit vérifiée “au mieux".
Ecrivons donc qu’on cherche ρk tel que ρk Hk ' I. En multipliant à droite par sk−1 = xk − xk−1 ,
on a
ρk Hk sk−1 ' sk−1
Mais Hk sk−1 ' yk−1 , car yk−1 = gk − gk−1 = ∇f (xk ) − ∇f (xk−1 ) ' Hk−1 (xk − xk−1 ) = Hk−1 sk−1
(cette dernière approximation résultant d’un développement limité à l’ordre 1) et Hk−1 ' Hk .
Ainsi, on obtient
ρk yk−1 ' sk−1 .
Ces considérations heuristiques nous amènent à définir ρk par
>
En minimisant ce trinôme ρ2 kyk−1 k2 − 2ρyk−1 sk−1 + ksk−1 k2 (on écrit que sa dérivée est nulle en
ρk ), on trouve
>
yk−1 sk−1
ρk =
kyk−1 k2
La méthode de Barzilai Borwein consiste à faire les itérations xk+1 = xk − ρk gk avec ce choix de
pas ρk . De manière surprenante, cette méthode s’avère être bien plus efficace que la méthode de
gradient à pas optimal ! A ce jour, on n’a toujours pas une explication complète de cette efficacité,
en dimension quelconque.
Symétriquement, on peut prendre aussi ρk = β1k où βk = argmin kyk−1 − βsk−1 k2 , ce qui donne
β>0
la variante
ksk−1 k2
ρk =
s>
k−1 yk−1
50
et notons que
d
f (x(t)) = −k∇f (x(t))k2
dt
autrement dit, f décroît le long de la trajectoire x(t) solution de (3.5). Dans les conditions de
convergence de la méthode de gradient ou de Newton, en fait, la trajectoire x(t) converge, lorsque
t → +∞, vers le minimiseur local x∗ . Cela est lié à la théorie de Lyapunov.
Cette observation constitue en fait la version continue des algorithmes de descente qu’on a vus
précédemment et donne un point de vue alternatif qui conduit à de nouvelles variantes d’algo-
rithmes.
– En effet, appliquons le schéma de discrétisation d’Euler explicite à (3.5) avec un pas ρk : on
obtient l’algorithme de descente
xk+1 = xk − ρk ∇f (xk )
– La méthode de Newton standard pour résoudre F (x) = ∇f (x) = 0 peut être vue de la
manière suivante : on écrit l’approximation à l’ordre 2
1
f (x + h) ' f (x) + ∇f (x)> h + h> Hf (x)h
2
et on minimise ce trinôme en h, ce qui donne h = Hf (x)−1 ∇f (x). Puis on remplace h par xk+1 −xk ,
ce qui conduit à
xk+1 = xk − Hf (xk )−1 ∇f (xk )
qui est bien la méthode de Newton.
– En appliquant le schéma de discrétisation d’Euler implicite à (3.5) avec un pas ρk , on obtient
xk+1 = xk − ρk ∇f (xk+1 )
qu’on appelle parfois la méthode d’Euler implicite linéarisée, ou encore modification de Levenberg-
Marquardt de la méthode de Newton. On note sur (3.6) que :
— si ρk est petit alors xk+1 ' xk − ρk ∇f (xk ) : méthode de descente ;
— si ρk est grand alors xk+1 ' xk − Hf (xk )−1 ∇f (xk ) : méthode de Newton. Autrement dit
(3.6) réalise un compromis entre la méthode de descente et la méthode de Newton.
La méthode d’Euler implicite linéarisée (3.6) est à comparer à la variante
qui s’appelle la “damped Newton method", car par rapport à la méthode standard de Newton on
ajoute le “damping" ρk . Cela permet de gérer les cas où la dérivée est quasi-singulière.
51
3.4.3 Méthode de gradient conjugué
L’algorithme du gradient conjugué ci-dessous est dû à Fletcher et Reeves (1964) :
— Initialisation : on choisit x0 ∈ IRn et d0 = −∇f (x0 ).
— Itération k (tant que ∇f (xk ) 6= 0) :
∇f (xk )> dk
ρk = −
d>
k Hf (xk )dk
xk+1 = xk + ρk dk
k∇f (xk+1 )k2
dk+1 = −∇f (xk+1 ) + dk
k∇f (xk )k2
— Critère d’arrêt.
L’idée générale de la méthode du gradient conjugué est de construire la suite (xk )k∈IN de façon à ce
que la suite (gk )k∈IN soit orthogonale, i.e., gk ⊥ gj pour k 6= j. En effet, grâce à cette propriété,
il existe un entier K 6 n − 1 tel que gK+1 = 0 (car les n + 1 vecteurs g0 , g1 , . . . , gn ne peuvent pas
être linéairement indépendants dans IRn ), i.e., AxK+1 = b et donc (par unicité) xK+1 = x̄. On va
voir que cet entier K se caractérise en analysant la suite des espaces de Krylov, définie ci-dessous.
Supposons que g0 = Ax0 −b 6= 0 (sinon, x0 = x̄ et il n’y a rien à faire). On pose u0 = g0 = Ax0 −b
puis on itère par A en posant
∀k ∈ IN∗ uk = Auk−1 = Ak u0
Notons que
∀k ∈ IN Uk ⊂ Uk+1 et AUk ⊂ Uk+1 . (3.7)
Soit
K = max{j ∈ IN | la famille (u0 , . . . , uj ) est libre}.
On a forcément K 6 n − 1. Par définition de K, uK+1 est combinaison linéaire des ui pour i 6 K.
Donc uK+2 = AuK+1 est combinaison linéaire des ui pour i 6 K + 1, donc combinaison linéaire
des ui pour i 6 K. Et ainsi de suite par récurrence. Par conséquent on a
52
autrement dit la suite d’espaces vectoriels Uk est strictement croissante jusqu’à k = K, puis
stationnaire.
Pour tout k ∈ IN∗ , soit xk ∈ Uk−1 l’unique minimiseur de f sur le sous-espace affine x0 + Uk−1 .
Comme tout point de Uk−1 est combinaison linéaire de u0 , . . . , uk−1 , on a
k−1
X
xk = argmin {f (x) | x ∈ x0 + Uk−1 } = argmin f x0 + a j uj a0 , . . . , ak−1 ∈ IR
j=0
Nous allons établir que, pour la fonction f quadratique considérée, la suite (xk )k∈IN converge en
exactement K + 1 itérations, ce qui est une propriété absolument remarquable ! Le problème est
de calculer de manière algorithmique les points xk . On va démontrer que ces points se calculent
itérativement par l’algorithme de Fletcher et Reeves donné en début de section. Pour cela, analysons
d’abord les propriétés de la suite (xk )k∈IN .
Faisons d’abord les remarques préliminaires suivantes. Comme xk ∈ x0 + Uk−1 , il peut s’écrire
k−1
X
xk = x0 + aj uj . Comme gk = Axk − b, on a donc
j=0
k−1
X k−1
X
gk = Axk − b = Ax0 − b + aj Auj = u0 + aj uj+1 .
| {z } |{z}
u0 j=0 j=0
uj+1
Par conséquent,
∀k ∈ IN gk ∈ Uk . (3.8)
Par ailleurs, pour tout k ∈ IN, on a gk+1 = Axk+1 − b = Ask + Axk − b, i.e.,
∀k ∈ IN sk ∈ Uk . (3.10)
Lemme 2. (i) La suite (gk )k∈IN est orthogonale, i.e., gk ⊥ gj pour k 6= j, ou de manière
équivalente,
∀k ∈ IN∗ gk ⊥ Vect(g0 , . . . , gk−1 ).
(ii) La suite (sk )k∈IN est A-orthogonale (on dit aussi A-conjuguée), i.e., Ask ⊥ sj pour k 6= j, ou
de manière équivalente,
Démonstration. Pour tout k ∈ IN∗ , par définition, xk minimise la fonction (a0 , . . . , ak−1 ) 7→
k−1
X
f x0 + aj uj sur IRk , donc la différentielle de cette fonction en xk (du moins, en le k-uplet de
j=0
coefficients de xk ) est nulle, ce qui donne
∀j ∈ {0, . . . , k − 1} gk> uj = 0
53
et donc
∀k ∈ IN∗ gk ⊥ Uk−1 . (3.11)
Montrons alors la propriété (i). D’après (3.8), on a gj ∈ Uj pour tout j ∈ IN. Soit k ∈ IN∗ .
Comme Uj ⊂ Uk−1 pour j 6 k − 1, on a donc Vect(g0 , . . . , gk−1 ) ⊂ Uk−1 . D’après (3.11), on a
gk ⊥ Uk−1 . On obtient donc (i).
Montrons maintenant (ii). D’après (3.11), on a gk+1 ⊥ Uk et gk ⊥ Uk−1 . Or, pour j 6 k − 1,
on a Uj ⊂ Uk−1 ⊂ Uk , donc gk+1 ⊥ Uj et gk ⊥ Uj . Or, d’après (3.9), on a Ask = gk+1 − gk , donc
Ask ⊥ Uj . Mais, d’après (3.10), sj ∈ Uj . D’où le résultat.
Montrons enfin (iii). Pour k = K, on a vu que UK+1 = UK . Mais alors, par (3.11), gK+1 ⊥ UK ,
et par (3.8), gK+1 ∈ UK+1 = UK , on en déduit donc que gK+1 = 0, i.e., AxK+1 = b et donc (par
unicité) xK+1 = x̄. C’est bien le premier entier pour lequel cette égalité arrive, car pour k 6 K la
famille (u0 , . . . , uk ) est libre et donc forcément gk 6= 0, i.e., xk 6= x̄.
D’après le lemme 2, (ii), on a, pour k 6= j, s>
k Asj = 0, or Asj = gj+1 − gj d’après (3.9), donc
s>
k (gj+1− gj ) = 0. Donc, en appliquant cette relation successivement pour j = k − 1, k − 2, . . . , 0,
on obtient
s> > >
k gk = sk gk−1 = · · · = sk g0 = αk
d’où
sk sk−1 gk
∀k ∈ {1, . . . , K} = + .
αk αk−1 kgk k2
En posant
kgk k2 sk αk
∀k ∈ {0, . . . , K} dk = − sk = , ρk = −
αk ρk kgk k2
on obtient bien :
• sk = ρk dk , i.e., xk+1 = xk + ρk dk ;
sk αk sk αk sk−1 αk gk αk ρk−1 dk−1 kgk k2
• dk = ρk = ρk αk = ρk αk−1 + ρk kgk k2 = ρk αk−1 − gk = −gk + kgk−1 k2 dk−1 ;
ce qui est l’algorithme écrit au début de la section, mais il reste encore à calculer ρk . Pour cela, on
note que gk+1 ⊥ gk d’après le lemme 2, (i), et gk+1 = gk + Ask d’après (3.9), donc 0 = gk> gk+1 =
kgk k2 + gk> Ask , et comme sk = ρk dk , on obtient
kgk k2
∀k ∈ {0, . . . , K} ρk = − .
gk> Adk
Or, d’après le lemme 2, (ii), la suite (dk )k∈IN est A-orthogonale (puisque la suite (sk )k∈IN l’est),
2
et comme gk = −dk + kgkgk−1 kk > > >
k2 dk−1 , on déduit que gk Adk = −dk Adk . Par ailleurs, gk dk =
2
− kgαkkk gk> sk = −kgk k2 car gk> sk = αk . D’où finalement
gk> dk
∀k ∈ {0, . . . , K} ρk = −
d>
k Adk
54
ce qui est la formule écrite en début de section. On peut noter qu’on a aussi
kgk k2
ρk =
d>k Adk
Conclusion. On a donc montré que, pour une fonction f quadratique, l’algorithme du gradient
conjugué de Fletcher et Reeves converge en exactement K + 1 6 n itérations.
On utilise cet algorithme pour des fonctions f non quadratiques, en pariant sur le fait que,
comme f est approchée à l’ordre deux par une fonction quadratique au voisinage d’un minimiseur
x̄, on s’attend à avoir un algorithme qui converge assez rapidement. Toutefois, à cause des termes
d’ordre supérieur, la convergence n’est plus exacte dans le cas général : il faut donc se donner un
critère d’arrêt, comme on l’a vu dans les sections précédentes.
Remarque 36. Notons que le pas ρk coïncide avec celui trouvé par la méthode du gradient à pas
optimal : en effet, en posant
ρk = argmin f (xk + ρdk )
ρ>0 | {z }
ϕ(ρ)
g > dk
on a 0 = ϕ0 (ρk ) = ∇f (xk + ρk dk )> dk ' gk> dk + ρk d>
k Hk dk , d’où ρk = − d> H
k
.
k k dk
Variante de Polak et Ribière (1969). Dans cette variante, on remplace la mise à jour de dk
par
∇f (xk+1 )> (∇f (xk+1 ) − ∇f (xk ))
dk+1 = −∇f (xk+1 ) + dk
k∇f (xk )k2
Lorsque f est quadratique, les deux algorithmes coïncident (car gk+1 ⊥ gk ). Toutefois, autant il
existe des résultats établissant la convergence de l’algorithme de Fletcher et Reeves, autant on ne
sait pas démontrer la convergence de l’algorithme de Polak et Ribière pour des larges classes de
fonctions !
L’algorithme du gradient conjugué de Polak et Ribière est le plus utilisé dans la pratique car
on constate qu’il converge et est plus performant pour de plus grandes classes de fonctions que
la version de Fletcher et Reeves (bien qu’on ne sache pas expliquer pourquoi, même 50 ans après
sa découverte !). On voit donc que, même sur des algorithmes classiques, on a encore des marges
d’amélioration et de découvertes potentielles...
Dans les deux algorithmes, pour être sûr d’avoir une direction de descente raisonnable, on peut
ajouter le test :
d>
k+1 ∇f (xk+1 )
6 −α < 0
kdk+1 kk∇f (xk+1 )k
où α > 0 est fixé, pas trop petit ; et si cela n’est pas vérifié, on prend simplement dk+1 =
−∇f (xk+1 ).
3.4.4 Conclusion
Dans cette section 3.4, on a vu diverses méthodes pour minimiser une fonction f : IRn → IR
suffisamment régulière. Les méthodes de type gradient consistent à suivre les lignes de gradient le
long des itérations : xk+1 = xk − ρk ∇f (xk ), où le pas ρk peut être astucieusement choisi. Dans les
méthodes de type Newton, l’itération s’écrit xk+1 = xk − Mk ∇f (xk ) où Mk est une matrice bien
55
choisie. Dans l’approche de gradient conjugué, la direction de descente est une combinaison linéaire
de ∇f (xk ) et de ∇f (xk−1 ), astucieusement choisie de façon à assurer des propriétés d’orthogonalité.
On peut combiner ces méthodes (comme on l’a vu dans les interprétations EDO) et en imaginer
d’autres, selon les problèmes considérés. Il ne faut pas hésiter à tester et à être créatif.
56
Chapitre 4
min f (x)
x∈C
min f (x)
x∈C
admet au moins un minimiseur (global), i.e., il existe au moins un point x∗ ∈ C tel que f (x∗ ) =
minx∈C f (x). On note aussi x∗ = argminx∈C f (x).
Lorsque E = IRn , et lorsque C est seulement fermé non vide (mais pas forcément borné), la
conclusion est encore vraie si de plus f est infinie à l’infini.
On n’a pas unicité du minimiseur en général (on parle donc de minimiseur local et de minimiseur
global). Pour avoir unicité, on peut supposer la convexité.
Théorème 24. Soit E un espace vectoriel et C ⊂ E un sous-ensemble convexe non vide. Si
f : E → IR est strictement convexe alors elle a au plus un minimiseur (global) sur C.
57
4.2 Conditions d’optimalité
4.2.1 Conditions d’optimalité du premier ordre sur un ensemble convexe
Théorème 25. Soit E un espace vectoriel normé, soit C ⊂ E un sous-ensemble convexe fermé
non vide, et soit f : E → IR une fonction Gateaux différentiable. Si x∗ ∈ C est un minimiseur
(local ou global) de f alors
∀x ∈ C df (x∗ ).(x − x∗ ) > 0
(on devrait plutôt écrire : f 0 (x∗ ; x − x∗ ) > 0 pour tout x ∈ C). Dans le cas où f est convexe, la
condition est nécessaire et suffisante (et x∗ est un minimiseur global).
Démonstration. Par convexité de C, x∗ + t(x − x∗ ) ∈ C pour tout t ∈ [0, 1], et pour t > 0 assez
petit on a f (x∗ + t(x − x∗ )) − f (x∗ ) > 0 car x∗ est un minimiseur au moins local. On divise par t
et on fait tendre t vers 0 pour obtenir la condition nécessaire.
Montrons qu’elle est suffisante si de plus f est convexe. Si f est convexe alors le graphe de f
est au-dessus de ses tangentes, donc f (x) > f (x∗ ) + df (x∗ ).(x − x∗ ), et donc, on en déduit que
f (x) > f (x∗ ), donc x∗ est un minimiseur (global).
Les conditions ci-dessus sont générales, sur un ensemble convexe, mais restent abstraites. Dans
la section suivante, on traite le cas E = IRn et un ensemble de contraintes général non convexe,
défini en termes d’égalités et d’inégalités de fonctions. Cela conduit à la notion de multiplicateurs
de Lagrange.
min f (x)
h(x)=0
p
X
λ∗0 ∇f (x∗ ) + λ∗i ∇hi (x∗ ) = 0
i=1
Si de plus les gradients ∇hi (x∗ ), i = 1, . . . , p sont indépendants (condition de qualification) alors
on peut de plus choisir λ∗0 = 1, et on a unicité des multiplicateurs de Lagrange.
58
Remarque 37. En supposant que λ est un vecteur colonne, et en identifiant dh(x∗ ) à la matrice
jacobienne, la condition des multiplicateurs de Lagrange s’écrit aussi sous la forme
Démonstration. La preuve est de nature géométrique et utilise le théorème des fonctions implicites.
On définit l’application “augmentée" F : IRn → IRp × IR par
Sur un dessin, représentons l’image de F , i.e., l’ensemble F (IRn ) : voir figure 4.1.
f (x)
F (IRn)
Ψ̃
F (x∗)
h(x) ∈ IRp
h=0
Im dF (x∗)
(droite verticale)
On veut minimiser f sur l’ensemble h = 0 qui, sur le dessin, correspond à l’axe des ordonnées :
on minimise donc f sur la “fibre" au-dessus de 0, qui intersecte l’ensemble F (IRn ). Le minimum est
obtenu au point le plus bas possible de cette fibre : c’est le point F (x∗ ) qui est noirci sur la figure.
Le point clé est de remarquer que, nécessairement, F (x∗ ) est au bord de l’ensemble F (IRn ) :
Ce point clé se trouve dans tout problème d’optimisation. En effet, être au bord, cela veut dire
qu’on ne peut pas faire mieux !
59
Notons deux choses : d’une part ce qu’on dit ici ne préjuge d’aucune considération spécifique
sur l’ensemble F (IRn ) (qu’on ne suppose pas fermé, ouvert, ou quoi que ce soit) ; tout ce qu’on
suppose c’est qu’un minimiseur existe, et alors, forcément, le point F (x∗ ) est au bord de l’ensemble
F (IRn ). D’autre part, la réciproque est fausse : être au bord de l’ensemble ne signifie pas qu’on a
un minimum ; on pourrait avoir un maximum, ou bien, rien du tout (ni minimum ni maximum).
N’oublions pas qu’on cherche ici une condition nécessaire d’optimalité.
Maintenant arrive l’étape du théorème des fonctions implicites : la fonction F est C 1 et, à
cause de la propriété de bord ci-dessus, la fonction F n’est pas localement surjective en x∗ (sinon,
il existerait une petite boule ouverte autour du point F (x∗ ) qui serait contenue dans F (IRn ) :
mais c’est faux ! cf figure de nouveau). Donc, par contraposée du théorème des fonctions implicites
(plus précisément, par contraposée du théorème de la submersion), la différentielle dF (x∗ ) : IRn →
IRp × IR n’est pas surjective, i.e.,
Im dF (x∗ ) ( IRp × IR
Comme on est en dimension finie, le sous-espace vectoriel strict Im dF (x∗ ) est donc contenu dans
un hyperplan 1 , et donc il existe un vecteur λ̃ ∈ IRp × IR \ {0} (non trivial ! sinon on ne dit rien...)
tel que λ̃∗ ⊥ Im dF (x∗ ), i.e., en notant λ̃∗ comme vecteur colonne,
Mais comme dF (x) = (df (x), dh(x)), en posant λ̃∗ = (λ∗ , λ∗0 ) avec λ∗ ∈ IRp et λ∗0 ∈ IR, on en
déduit la relation des multiplicateurs de Lagrange.
Sous la condition de qualification, on a forcément λ∗0 6= 0, car sinon, si λ∗0 = 0 on aurait une
relation de dépendance linéaire entre les gradients ∇hi (x∗ ). On peut alors normaliser le (p+1)-uplet
des multiplicateurs de Lagrange de sorte que λ∗0 = 1, grâce à la remarque 38 ci-dessous. De plus
on a unicité des multiplicateurs de Lagrange (λ∗1 , . . . , λ∗p ) : en effet si un p-uplet (λ̂∗1 , . . . , λ̂∗p ) vérifie
Pp Pp
∇f (x∗ ) + i=1 λ̂∗i ∇hi (x∗ ) = 0, alors en retranchant la relation ∇f (x∗ ) + i=1 λ∗i ∇hi (x∗ ) = 0, on
Pp
obtient i=1 (λ̂∗i − λ∗i )∇hi (x∗ ) = 0, et par indépendance des gradients on obtient λ̂∗i = λ∗i pour
tout i.
Remarque 38. Il est important de noter que le multiplicateur de Lagrange λ̃∗ = (λ∗ , λ∗0 ) construit
ci-dessus, d’une part, est non trivial (sinon on écrit 0 = 0 ! donc rien...), et d’autre part, est défini
à scalaire multiplicatif près, autrement dit, pour tout α 6= 0, αλ̃∗ est aussi un multiplicateur de
Lagrange.
Le réel λ∗0 est appelé multiplicateur de Lagrange associé au coût.
Pour i ∈ {1, . . . , p}, le réel λ∗i est appelé multiplicateur de Lagrange associé à la contrainte
hi = 0.
On a deux cas possibles :
• Si λ∗0 6= 0, quitte à multiplier le multiplicateur de Lagrange λ̃∗ par 1/λ∗0 , on peut supposer
que λ∗0 = 1. Ce cas s’appelle le cas normal.
• Mais il se peut que λ∗0 = 0 : on appelle ce cas le cas anormal.
Bien entendu, le cas anormal n’arrive pas sous la condition de qualification. Mais sinon, il pourrait
arriver : par exemple considérons le problème d’optimisation (certes trivial) avec n = 1, f (x) = x
et h(x) = x2 ; alors on a un cas anormal. On comprend ici que cela vient du fait que l’ensemble
h = 0 est un point isolé. Dans le cas général, cela arrive lorsque l’ensemble h = 0 n’est pas une
sous-variété de IRn au point x∗ (point singulier).
1. Notons que ce fait pourrait échouer en dimension infinie : on peut avoir un sous-espace vectoriel strict qui soit
partout dense ! En dimension infinie, il faut alors trouver des hypothèses qui impliquent que Im dF (x∗ ) est de plus
fermé ; auquel cas l’argument de séparation marche encore...
60
Remarque 39. La condition de qualification s’exprime de manière équivalente en disant que la
différentielle dh(x∗ ) : IRn → IRp est surjective.
Sous cette condition, le théorème des fonctions implicites implique que l’ensemble des contraintes
h = 0 est une sous-variété de IRn (de codimension p).
C’est le sens des conditions de qualification : lorsque l’ensemble de contraintes est une sous-
variété, alors il existe un multiplicateur de Lagrange normal ; mais lorsque l’ensemble de contraintes
a un point singulier, on peut avoir un multiplicateur de Lagrange anormal.
Il faut prendre garde et ne pas oublier les multiplicateurs anormaux, lorsqu’on recherche les
solutions optimales !
ce qui donne, vu que ∇f (0, 0) = (1 1)> , ∇h1 (0, 0) = (−2 0)> , ∇h2 (0, 0) = (−4 0)> : λ∗0 = 0 et
λ∗1 = −2λ∗2 . Comme on s’y attend, il reste un degré de liberté, vu que le multiplicateur de Lagrange
est défini à scalaire multiplicatif près. On peut prendre par exemple λ∗0 = 0, λ∗1 = −2, λ∗2 = 1. Quoi
qu’il en soit, c’est un multiplicateur anormal.
Remarque 40. Lorsque toutes les fonctions hi , i ∈ {1, . . . , p} sont affines, il existe toujours au
moins un multiplicateur normal. En effet, dans ce cas, comme on suppose qu’il existe un minimiseur
x∗ , si les contraintes hi ne sont pas indépendantes alors certaines d’entre elles sont redondantes
(d’où l’existence d’un multiplicateur anormal) ; mais sans perte de généralité on peut ignorer les
contraintes redondantes et se ramener au cas qualifié, d’où l’existence d’un multiplicateur normal.
Mais dans ce cas où les contraintes ne sont pas indépendantes, on n’a évidemment pas unicité du
multiplicateur.
∂L ∗ ∗ ∗
Avec cette fonction, la condition de multiplicateurs de Lagrange s’écrit ∂x (x , λ0 , λ ) = 0. Par
ailleurs, on a h(x∗ ) = 0, mais h = ∂L
∂λ . Donc finalement :
∂L ∗ ∗ ∗ ∂L ∗ ∗ ∗
(x , λ0 , λ ) = 0 et (x , λ0 , λ ) = 0
∂x ∂λ
autrement dit, (x∗ , λ∗ ) est un point extrémal de la fonction (x, λ) 7→ L(x, λ∗0 , λ), ce qui veut dire
aussi :
• x∗ est un point extrémal de la fonction x 7→ L(x, λ∗0 , λ∗ ) ;
• λ∗ est un point extrémal de la fonction λ 7→ L(x∗ , λ∗0 , λ).
61
Remarque 41. Si on suppose de plus que f est convexe et que les fonctions hi sont affines, alors
∗ ∗ ∗ ∗
L est convexe par rapport à x, donc la condition ∂L
∂x (x , λ0 , λ ) = 0 est équivalente au fait que x
∗ ∗
minimise la fonction L(·, λ0 , λ ) :
En général, toutefois, il n’est pas vrai que x∗ minimise cette fonction. On reviendra toutefois sur
ce point plus tard, d’une part dans la section [Link] sur la méthode de Lagrangien augmenté et
d’autre part dans la section 4.3.2 sur les méthodes duales.
On peut aussi d’ores-et-déjà noter que, sans aucune hypothèse supplémentaire de convexité,
λ∗ est un maximiseur de la fonction λ 7→ L(x∗ , λ∗0 , λ), i.e., L(x∗ , λ∗0 , λ) 6 L(x∗ , λ∗0 , λ∗ ) pour tout
λ ∈ IRp . On a même beaucoup mieux que cela : en fait on a l’égalité L(x∗ , λ∗0 , λ) = L(x∗ , λ∗0 , λ∗ ) =
λ∗0 f (x∗ ) car h(x∗ ) = 0. Mais on verra plus loin que la propriété générale, lorsqu’il y a des contraintes
d’inégalité, est la propriété de maximisation.
Lorsqu’on a un multiplicateur de Lagrange normal (i.e., λ0 > 0), qu’on normalise à λ0 = 1,
souvent on note
p
X
L(x, λ) = L(x, 1, λ) = f (x) + λi hi (x).
i=1
min f (x)
h(x)=0
g(x)60
où f : IRn → IR, g : IRn → IRq et h : IRn → IRp sont C 1 . Ici, la notation g(x) 6 0 signifie que
gj (x) 6 0 pour tout j = 1, . . . , q.
p
X q
X
λ∗0 ∇f (x∗ ) + λ∗i ∇hi (x∗ ) + µ∗j ∇gj (x∗ ) = 0
i=1 j=1
∗ ∗
h(x ) = 0, g(x ) 6 0
λ∗0 > 0, µ∗j > 0 ∀j ∈ {1, . . . , q}
µ∗j gj (x∗ ) =0 ∀j ∈ {1, . . . , q} (conditions de complémentarité)
Si de plus les gradients ∇hi (x∗ ), i ∈ {1, . . . , p}, ∇gj (x∗ ), j ∈ I(x∗ ) (indices des contraintes
actives) sont tous linéairement indépendants (famille libre), alors on peut de plus supposer que λ0 =
1 ci-dessus, et on a unicité des multiplicateurs de Lagrange. On parle de condition de qualification.
La contrainte gj 6 0 est dit active (ou saturée) si gj (x∗ ) = 0, et inactive (ou non saturée)
si gj (x∗ ) < 0. On note I(x∗ ) l’ensemble des indices actifs, i.e., des indices correspondant à une
contrainte active.
Les conditions de complémentarité sont une manière équivalente de dire que, pour j ∈ {1, . . . , q},
si la contrainte gj (x∗ ) est inactive en x∗ , i.e., si gj (x∗ ) < 0, alors µ∗j = 0 : le multiplicateur de
62
Lagrange correspondant est nul. Au contraire lorsque la contrainte est active, i.e., si gj (x∗ ) = 0,
alors la relation µ∗j gj (x∗ ) = 0 est bien vérifiée (on pourrait avoir malgré tout µ∗j = 0, mais en
général on aura µ∗j > 0).
Donnons deux démonstrations du théorème KKT.
Première démonstration de KKT. Le problème ci-dessus est équivalent à un problème d’optimi-
sation comportant uniquement des contraintes d’égalité : ce sont les contraintes hi (x) = 0 d’une
part, avec i = 1, . . . , q, et d’autre part les contraintes gj (x) = 0 pour les indices actifs j ∈ I(x∗ ).
Bien sûr, comme on ne connaît pas à l’avance le minimiseur x∗ , on ne sait pas non plus à
l’avance quelles sont les contraintes actives ! Tout cela est théorique.
Mais il n’empêche que, de fait, x∗ est aussi une solution optimale du problème avec contraintes
d’égalité
min f (x) (4.1)
H(x)=0
∗
où H : IRn → IRp × IRcard(I(x ))
est donné par
h(x)
H(x) = . (4.2)
gj (x), j ∈ I(x∗ )
Notons que le problème est qualifié si et seulement si l’application linéaire dH(x∗ ) : IRn → IRp ×
∗
IRcard(I(x )) est surjective (condition sous laquelle le théorème des fonctions implicites implique que
l’ensemble {x ∈ IRn | H(x) = 0} est une sous-variété de IRn ).
Par le théorème des multiplicateurs de Lagrange, il existe donc λ∗0 ∈ IR et un vecteur colonne
∗
Ψ ∈ IRp × IRcard(I(x )) , de coordonnées successives λ∗1 , . . . , λ∗p , (µ∗j )j∈I(x∗ ) , tels que
La preuve ci-dessus est une preuve par “argument de séparation", en conservant des considéra-
tions de signe.
63
f (x)
F (IRn)
Ψ̃
F (x∗)
g(x) ∈ IRq
g 60
Im dF (x∗)
(espace tangent à F (IRn) en F (x∗))
Figure 4.2 – Image de l’application augmentée F (x) = (h(x), g(x), f (x)) (on ne représente pas les
coordonnées h), et visualisation de la solution optimale : F (x∗ ) est au bord de l’ensemble F (IRn ).
Notons bien qu’on aurait pu prendre la convention opposée, et choisir λ∗0 6 0 : dans ce cas on
aurait aussi µ∗j 6 0.
Il existe d’autres preuves, comme celle ci-dessous qui est intéressante car elle ouvre la voie aux
techniques de pénalisation (qu’on verra plus loin) :
Deuxième démonstration de KKT. Soit R > 0 quelconque, on note B̄(x∗ , R) la boule fermée de
IRn de centre x∗ et de rayon R. Pour tout ε > 0, on définit la “fonction pénalisée" :
p q
1 X 1 X
fε (x) = f (x) + hi (x)2 + max(gj (x), 0)2 + kx − x∗ k2 .
2ε i=1 2ε j=1
1
L’idée ici est que, si on prend x tel que hi (x) 6= 0 alors le terme 2ε hi (x)2 devient très grand lorsque
1
ε est petit. De même, si x est tel que gj (x) > 0 alors le terme 2ε max(gj (x), 0)2 devient très grand
lorsque ε est petit. On s’attend donc à ce que, lorsque ε → 0, "le" (ou, "un") minimiseur xε de
fε tend vers un point vérifiant les contraintes. Le terme supplémentaire kx − x∗ k2 force ce point à
être égal au point désiré x∗ .
Avec cette idée, faisons alors rigoureusement la preuve. On considère le problème
min fε (x).
x∈B̄(x∗ ,R)
64
Tout d’abord, ce problème admet au moins un minimiseur xε ∈ B̄(x∗ , R), car fε est continue sur
le compact B̄(x∗ , R).
Comme xε ∈ B̄(x∗ , R), à sous-suite près on peut supposer que xε → x̄ ∈ B̄(x∗ , R) lorsque
ε → 0. On va démontrer que, en fait, x̄ = x∗ .
Tout d’abord, comme xε minimise fε , on a fε (xε ) 6 fε (x∗ ), et par ailleurs, comme h(x∗ ) = 0
et g(x∗ ) 6 0, on observe que fε (x∗ ) = f (x∗ ). Ainsi, on a fε (xε ) 6 f (x∗ ), et donc
p
X q
X
hi (xε )2 + max(gj (xε ), 0)2 6 2ε f (x∗ ) − f (xε ) − kxε − x∗ k2 6 2εM
i=1 j=1
pour un M > 0 car l’expression entre parenthèses est bornée (on est sur un compact). Donc, lorsque
ε → 0, on a hi (xε ) → 0 et max(gj (xε ), 0) → 0. En passant à la limite, comme xε → x̄ on obtient
donc
hi (x̄) = 0 ∀i ∈ {1, . . . , p} gj (x̄) 6 0 ∀j ∈ {1, . . . , q}
autrement dit, x̄ vérifie les contraintes.
Par ailleurs, par définition de fε on a f (xε ) + kxε − x∗ k2 6 fε (xε ) (puisqu’on ajoute des carrés),
et on a vu que fε (xε ) 6 f (x∗ ), donc f (xε ) + kxε − x∗ k2 6 f (x∗ ), et en passant à la limite on
obtient f (x̄) + kx̄ − x∗ k2 6 f (x∗ ). Mais comme x∗ est un minimiseur du problème contraint et que
x̄ est un point qui vérifie les contraintes, on doit forcément avoir f (x∗ ) 6 f (x̄), ce qui donne donc
finalement f (x̄) + kx̄ − x∗ k2 6 f (x̄) et donc x̄ = x∗ .
On note que, comme ce raisonnement a été fait pour toute sous-suite convergente de xε , fina-
lement, xε → x∗ (pas seulement à sous-suite près).
Ecrivons maintenant les conditions nécessaires d’optimalité pour fε . Comme xε minimise fε sur
la boule B̄(x∗ , R), et comme xε est dans l’intérieur de la boule lorsque ε est assez petit (puisque
xε → x∗ ), la condition nécessaire d’optimalité est
∇fε (xε ) = 0
d
(notons ici qu’on a utilisé la formule dx max(x, 0)2 = 2 max(x, 0), facile à montrer)
Cette égalité commence à ressembler à la relation de multiplicateurs de Lagrange ! En effet si
on pose
hi (xε ) max(gj (xε ), 0)
λεi = , µεj =
ε ε
on a
Xp Xq
∇f (xε ) + λεi ∇hi (xε ) + µεj ∇gj (xε ) + 2(xε − x∗ ) = 0.
i=1 j=1
On veut maintenant faire tendre ε vers 0 dans cette égalité. Pour faire cela, on pose
65
ε
ψ
La famille de vecteurs kψ ε k est de norme 1, donc à sous-suite près, elle converge lorsque ε → 0
vers un vecteur (λ0 , λ , µ ) ∈ IR × IRp × IRq de norme 1. En notant que kψ ε k > 1 et en passant à
∗ ∗ ∗
la limite on a donc
p
X q
X
λ∗0 ∇f (x∗ ) + λ∗i ∇hi (x∗ ) + µ∗j ∇gj (x∗ ) = 0
i=1 j=1
ce qui est la relation des multiplicateurs de Lagrange, et de plus par la construction ci-dessus on
a bien obtenu µ∗j > 0, ce qui était la chose nouvelle à obtenir.
Notons que si gj (x∗ ) < 0 (contrainte inactive) alors gj (xε ) < 0 pour ε assez petit et donc µεj = 0
et donc µ∗j = 0 par passage à la limite.
Comme dans la preuve des multiplicateurs de Lagrange, si les gradients sont indépendants alors
λ∗0 6= 0 (par l’absurde).
Formulation Lagrangienne. On peut formuler les conditions KKT sous la forme suivante. On
définit la fonction L : IRn × IR × IRp × IRq+ (Lagrangien) par
p
X q
X
L(x, λ0 , λ, µ) = λ0 f (x) + λi hi (x) + µj gj (x).
i=1 j=1
∂L ∗ ∗ ∗ ∗
Avec cette fonction, la condition KKT s’écrit ∂x (x , λ0 , λ , µ ) = 0. Par ailleurs, on a d’une part
h(x∗ ) = 0, mais h = ∂L
∂λ , et d’autre part on a
∂L ∗ ∗ ∗ ∗ ∗
∂µj (x , λ0 , λ , µ ) = gj (x ) 6 0 et qui est nul pour
j ∈ I(x∗ ). Finalement :
∂L ∗ ∗ ∗ ∗ ∂L ∗ ∗ ∗ ∗ ∂L ∗ ∗ ∗ ∗ 60
(x , λ0 , λ , µ ) = 0, (x , λ0 , λ , µ ) = 0, (x , λ0 , λ , µ )
∂x ∂λ ∂µj =0 ∀j ∈ I(x∗ )
autrement dit, (x∗ , λ∗ , µ∗ ) est un point extrémal de la fonction (x, λ, µ) 7→ L(x, λ∗0 , λ, µ) sur IRn ×
IRp × IRq+ , ce qui veut dire aussi :
• x∗ est un point extrémal de la fonction x 7→ L(x, λ∗0 , λ∗ , µ∗ ) sur IRn ;
• (λ∗ , µ∗ ) est un point extrémal de la fonction (λ, µ) 7→ L(x∗ , λ∗0 , λ, µ) sur IRp × IRq+ .
Remarque 43. En fait, comme évoqué dans la remarque 41, (λ∗ , µ∗ ) est un maximiseur de la
fonction (λ, µ) 7→ L(x∗ , λ∗0 , λ, µ). En effet, on a h(x∗ ) = 0, et comme g(x∗ ) 6 0, on a aussi
µ> g(x∗ ) 6 0 pour tout µ ∈ IRq+ , et donc
L(x∗ , λ∗0 , λ, µ) = f (x∗ ) + λ> h(x∗ ) + µ> g(x∗ ) 6 f (x∗ ) = L(x∗ , λ∗ , µ∗ ) ∀(λ, µ) ∈ IRp × IRq+
(pour la dernière égalité L(x∗ , λ∗ , µ∗ ) = f (x∗ ), on a utilisé le fait que h(x∗ ) = 0 et l’égalité
(µ∗ )> g(x∗ ) = 0 qui vient des conditions de complémentarité), ce qui montre la propriété de maxi-
misation du Lagrangien par rapport à (λ, µ) ∈ IRp × IRq+ .
Par ailleurs, et comme évoqué aussi dans la remarque 41, x∗ n’est pas forcément un minimiseur
de la fonction x 7→ L(x, λ∗0 , λ∗ , µ∗ ), mais il l’est sous des conditions de convexité (ce qui ouvre
aux méthodes dites de dualité qu’on va voir plus loin), comme le montre la preuve du résultat
ci-dessous qui, du coup, donne une réciproque aux conditions KKT. De plus, dans ce cas, on n’a
pas besoin de condition de qualification.
66
Proposition 4. Par rapport au théorème KKT, on suppose de plus que f est convexe, que les
fonctions gj sont convexes et les fonctions hi sont affines. Alors, x∗ est un minimiseur si et
seulement si on a les conditions KKT avec un multiplicateur normal (λ∗0 = 1).
Démonstration. Avec les hypothèses de convexité, L est convexe par rapport à x (notons que
∗ ∗ ∗
µj > 0), donc, comme dans la remarque 41, la condition ∂L
∂x (x , λ , µ ) = 0 est équivalente au fait
que x∗ minimise la fonction L(·, λ∗ , µ∗ ) :
L(x∗ , λ∗ , µ∗ ) 6 L(x, λ∗ , µ∗ ) ∀x ∈ IRn . (4.3)
Soit alors x ∈ IRn un point vérifiant les contraintes h(x) = 0 et g(x) 6 0. Comme µ∗j > 0, on a
Pp Pq
donc i=1 λ∗i hi (x) + j=1 µ∗j gj (x) 6 0, ce qui implique que L(x, λ∗ , µ∗ ) 6 f (x). Ainsi, en notant
que L(x∗ , λ∗ , µ∗ ) = f (x∗ ) (puisque h(x∗ ) = 0 et µ∗j gj (x∗ ) = 0), l’inégalité (4.3) donne
f (x∗ ) 6 f (x) ∀x ∈ IRn ,
ce qui montre que x∗ est minimiseur du problème avec contraintes.
Remarque 44 (Contraintes redondantes et multiplicateur anormal.). Quand on résout un pro-
blème d’optimisation sous contraintes d’égalité et d’inégalité, prenons garde à éviter les contraintes
redondantes, qui créent un multiplicateur anormal (i.e., tel que λ∗0 = 0).
Cela signifie que, en notant C = {x ∈ IRn | h(x) = 0, g(x) 6 0}, on essaie d’éliminer, autant
que possible, toutes les contraintes qui sont redondantes : par exemple, si le problème comporte
les contraintes x1 + x2 = 1 et 2x1 + 2x2 = 2, on ne garde bien sûr que la première.
Si des contraintes sont redondantes, il existe toujours un multiplicateur de Lagrange anormal !
Par exemple, considérons le problème de minimiser f : IRn → IR (de classe C 1 , quelconque) sous les
contraintes d’égalité h1 (x) = b> x−1 = 0 et h2 (x) = 2b> x−2 = 0 (où b ∈ IRn \{0} est fixé), qui sont
évidemment redondantes. On a ∇h1 = b et ∇h2 = 2b = 2∇h1 , donc 0 × ∇f (x) + 2∇h1 − ∇h2 = 0
et donc l’égalité des multiplicateurs de Lagrange est vérifiée avec le multiplicateur anormal (λ∗0 =
0, λ∗1 = 2, λ∗2 = 1).
(LI) Les gradients ∇hi (x∗ ), i ∈ {1, . . . , p}, ∇gj (x∗ ), j ∈ I(x∗ ) (indices des contraintes actives)
sont tous linéairement indépendants (famille libre).
Ici, le sigle (LI) signifie “linéairement indépendant". Comme on l’a vu dans le théorème 27, la
condition (LI) implique non seulement qu’on est dans le cas normal (i.e., on peut choisir λ0 = 1
par renormalisation), mais implique aussi l’unicité du multiplicateur de Lagrange (qui de plus est
normal).
Dans la démonstration de KKT vue ci-dessus, on voit qu’on aurait pu affaiblir cette hypothèse
et toutefois maintenir la conclusion de multiplicateur normal (mais pas forcément l’unicité !). Pour
démontrer par l’absurde que λ∗0 6= 0, il suffit de faire l’hypothèse plus faible suivante :
(PI) Les gradients ∇hi (x∗ ), i ∈ {1, . . . , p}, ∇gj (x∗ ), j ∈ I(x∗ ) (indices des contraintes actives)
card(I(x∗ ))
sont IRp × IR+ -indépendants, au sens suivant :
card(I(x∗ ))
∀λ = (λ1 , . . . , λp ) ∈ IRp ∀µ = (µj )j∈I(x∗ ) ∈ IR+
p
X X
λi ∇hi (x∗ ) + µj ∇gj (x∗ ) = 0 =⇒ λ = 0, µ = 0.
i=1 j∈I(x∗ )
67
Ici, le sigle (PI) signifie, avec un léger abus, “positivement indépendant".
La condition (PI) est suffisante, dans le théorème KKT, pour assurer l’existence d’un multipli-
cateur normal (λ∗0 = 1). Par contre, elle ne garantit pas l’unicité du multiplicateur de Lagrange.
Elle est plus générale que la condition (LI), i.e., (LI) ⇒ (PI). Par exemple si p = 0 (i.e., si on n’a
que des contraintes d’inégalité), on peut avoir une infinité de vecteurs qui sont “positivement indé-
pendants" au sens ci-dessus : il suffit de prendre des vecteurs qui pointent tous dans le quadrant
positif (c’est-à-dire, dont toutes les coordonnées sont positives) ; et bien sûr, dès que ces vecteurs
sont trop nombreux ils ne peuvent pas être linéairement indépendants.
Il existe de nombreuses autres conditions dans la littérature, assurant l’existence d’un multi-
plicateur normal. Ci-dessous, on en donne deux, qui sont très connues (toutefois, moins générales
que (PI)) : la condition de Mangasarian-Fromovitz, et la condition de Slater.
(MF) Les gradients ∇hi (x∗ ), i ∈ {1, . . . , p}, sont linéairement indépendants, et il existe d ∈
IRn \ {0} tel que
h∇gj (x∗ ), di < 0, on en déduit que µj = 0 pour tout j ∈ I(x∗ ). Puis on déduit que λi = 0 pour
tout i ∈ {1, . . . , p} puisque les vecteurs ∇hi (x∗ ) sont linéairement indépendants.
Condition de Slater. Cette condition, qui a l’avantage d’être facile à vérifier, concerne plus
spécifiquement les problèmes d’optimisation avec contraintes d’égalité affines et contraintes d’in-
égalité convexes :
(S) Les fonctions hi , i ∈ {1, . . . , p}, sont affines, les fonctions gj , j ∈ {1, . . . , q}, sont convexes,
et il existe x ∈ IRn tel que h(x) = 0 et g(x) < 0.
68
Remarque 45. On peut légèrement généraliser la condition de Slater de la manière suivante, en
distinguant, parmi les fonctions gj , celles qui sont affines de celles qui sont convexes non affines. La
condition est alors : les fonctions hi , i ∈ {1, . . . , p}, sont affines et les fonctions gj , j ∈ {1, . . . , q},
sont convexes ; il existe x ∈ IRn tel que h(x) = 0 et gj (x) < 0 si gj n’est pas affine, pour j = 1, . . . , q.
Exemple 9. Considérons le problème
1
min kx − y0 k2 | x1 > 0, x2 > 0, x1 + x2 6 1
2
avec y0 = (1, 1/2). Bien sûr, la solution est évidente (faire un dessin) : la solution x∗ est le projeté
orthogonal de y0 sur le triangle hachuré. On retrouve toutefois par le calcul, après résolution des
conditions KKT :
3 1 1
x∗1 = , x∗2 = , λ∗0 = 1, µ∗1 = µ∗2 = 0, µ∗3 = .
4 4 4
i.e., le problème de minimiser la fonction f (x) = 12 x> Ax − b> x + c sous les k contraintes d’égalité
Cx = d : en notant L1 , . . . , Lk ∈ M1,n (IR) les k lignes de la matrice C, ces k contraintes affines
sont
hi (x) = Li x − di = 0, i = 1, . . . , k.
On suppose C surjective, i.e., rg(C) = k, ce qui est logique car on souhaite que l’ensemble des
x tels que Cx = d soit non trivial ! Cela veut exactement dire que les k contraintes d’égalité
sont indépendantes, donc la condition de qualification est vérifiée : les vecteurs ∇hi (x) = L> i ,
i = 1, . . . , k, sont linéairement indépendants.
Le problème (4.4) a une unique solution x car f est strictement convexe et infinie à l’infini, et
l’ensemble des contraintes est convexe. D’après la proposition 4, l’unique solution x du problème est
caractérisée par la condition de multiplicateurs de Lagrange (qui est alors une condition nécessaire
et suffisante) : il existe λ = (λ1 , . . . , λk ) ∈ IRk tel que
k
X
∇f (x) + λi ∇hi (x) = 0
i=1
c’est-à-dire
Ax + C > λ = b
Cx = d
ce qui s’écrit sous forme matricielle
C>
A x b
= (4.5)
C 0 λ d
| {z }
M
69
Notons bien que (4.5) est une condition nécessaire et suffisante pour que x soit solution du problème
(4.4), avec λ comme multiplicateur de Lagrange associé (et λ0 = 1 car le problème est qualifié).
Notons aussi que la matrice M ∈ Mn+k(IR) apparaissant
dans (4.5) est (symétrique) inversible.
n k y 0
En effet, soit (y, µ) ∈ IR × IR tel que M = , alors on a Ay + C > µ = 0 et Cy = 0, donc
µ 0
y
0 = y> µ> M = y > Ay + 2µ> Cy = y > Ay
µ
et donc y = 0 puisque A est symétrique définie positive. On en déduit que C > µ = 0, et comme
C > est injective (car C est surjective), on a aussi µ = 0. Donc M est inversible.
En particulier, le multiplicateur de Lagrange normal λ associé à x est unique.
Cas particulier : moindres carrés contraints. Dans une séance précédente, nous avons vu
comment résoudre le problème des moindres carrés
avec A ∈ Mm,n (IR), b ∈ IRn , qui est un problème dont l’importance est cruciale dans un grand
nombre d’applications. On a vu que ce problème admet une unique solution lorsque A est injective
(ce qui impose m > n : la matrice A a plus de lignes que de colonnes), donnée par la pseudo-inverse
On s’intéresse ici à ajouter des contraintes d’égalité affines dans cet important problème :
où C ∈ Mk,n (IR) est surjective et d ∈ IRk , i.e., on ajoute k contraintes d’égalité affines indé-
pendantes. Bien entendu, on peut mettre un 21 devant la norme dans (4.6) : cela ne change rien
au problème de minimisation, et alors la fonctionnelle f à minimiser est du type précédent :
f (x) = 12 x> A> Ax − b> Ax + 12 kbk2 (la matrice A> A est symétrique définie positive).
Avec ce qu’on a vu ci-dessus, il existe une unique solution optimale x ∈ IRn , qui admet un
unique multiplicateur de Lagrange normal λ ∈ IRk , et le couple (x, λ) est l’unique solution du
système
>
A A C>
>
x A b
=
C 0 λ d
Application : distance d’un point à un sous-espace affine. Dans le problème (4.6), prenons
m = n et A = In . On a alors le problème
min kx − bk2
x∈IRn
Cx=d
qui est exactement le problème de trouver, pour un point b ∈ IRn donné, quel est le point x du sous-
espace affine F = {x ∈ IRn | Cx = d} qui minimise la distance à b. Bien sûr, géométriquement,
ce point x est exactement le projeté orthogonal du point b sur F : x = PF (b).
70
D’après ce qu’on a vu ci-dessus, cet unique point est caractérisé par les équations
x + C > λ = b, Cx = d
On peut noter que, comme C est surjective, sa pseudo-inverse est C # = C > (CC > )−1 .
Notons que la projection PF ci-dessus est une application affine. Si d = 0, la formule ci-dessus
donne x = In − C > (CC > )−1 C b qui est exactement la projection orthogonale de b sur le sous-
espace F = ker(C). On a donc obtenu, en prime, le résultat suivant : le projecteur orthogonal sur
le sous-espace vectoriel ker(C) est
min kxk2 .
Ax=b
d’après ce qu’on a dit précédemment, il existe une unique solution, caractérisée par les équations
x + A> λ = 0 et Ax = b, ce qui conduit à
Notons que, comme A est surjective, sa pseudo-inverse est A# = A> (AA> )−1 .
Comme expliqué dans le chapitre sur les moindres carrés, toutefois, ce cas est moins intéressant.
Le cas intéressant en pratique est de vouloir résoudre Ax = b avec A injective (A ayant plus de
lignes que de colonnes), et comme ce système est sans solution on résout le problème des moindres
carrés pour trouver la meilleure “solution approchée" possible.
71
Considérons le problème d’optimisation
min f (x)
h(x)=0
g(x)60
où f : IRn → IR, g : IRn → IRq et h : IRn → IRp sont C 2 . On a vu précédemment qu’un tel problème
peut toujours se ramener au problème (4.1) avec contraintes d’égalité seulement, en introduisant
les indices actifs et l’application H définie par (4.2). Ainsi, on se ramène au problème
min f (x),
H(x)=0
Soit (x∗ , λ∗0 , λ∗ , µ∗ ) ∈ IRn × IR+ × IRp × IRq+ une solution extrémale de KKT, i.e.,
∂L ∗ ∗ ∗ ∗ ∂L ∗ ∗ ∗ ∗ ∂L ∗ ∗ ∗ ∗
(x , λ0 , λ , µ ) = 0, (x , λ0 , λ , µ ) = 0, (x , λ0 , λ , µ ) = 0 ∀j ∈ I(x∗ ),
∂x ∂λ ∂µj
autrement dit, vérifiant la condition nécessaire d’optimalité du premier ordre.
On veut maintenant donner des conditions de second ordre.
On suppose que ce problème est qualifié, au sens du théorème 26, c’est-à-dire que les gradients
des composantes de H en x∗ , qui sont l’ensemble des vecteurs ∇hi (x∗ ) pour i = 1, . . . , p et ∇gj (x∗ )
pour j ∈ I(x∗ ), sont linéairement indépendants ; de manière équivalente, dH(x∗ ) est surjective
(cette condition pourrait être affaiblie). L’ensemble H = 0 est alors une sous-variété localement
autour de x∗ , et son espace tangent est
p
!
\ \ \
ker dH(x∗ ) = {d ∈ IRn | dH(x∗ ).d = 0} = ker dhi (x∗ ) ker dgj (x∗ ) .
i=1 j∈I(x∗ )
Dans ces conditions, on est dans le cas normal, donc quitte à renormaliser on peut toujours supposer
que λ∗0 = 1, et de plus il existe un unique multiplicateur de Lagrange (λ∗ , µ∗ ) ∈ IRp × IRq+ . On a
alors le théorème suivant.
Théorème 28. Sous les hypothèses ci-dessus :
— Condition nécessaire d’optimalité du second ordre : si x∗ est un minimiseur, alors
∂2L ∗ ∗ ∗
(x , λ , µ ).(d, d) > 0 ∀d ∈ ker dH(x∗ ).
∂x2
— Condition suffisante d’optimalité locale du second ordre : si
∂2L ∗
(x , λ, µ).(d, d) > 0 ∀d ∈ ker dH(x∗ ) \ {0}
∂x2
alors x∗ est un minimiseur local.
Démonstration. On peut admettre cette preuve dans un premier temps, car elle est plus difficile. La
preuve fait en effet appel à un concept de géométrie différentielle un peu élaboré : le concept de déri-
vée seconde intrinsèque. Pour l’expliquer, considérons l’application augmentée F (x) = (f (x), H(x)).
On a vu que les équations KKT sont équivalentes à la relation Ψ̃ ⊥ Im dF (x∗ ) (i.e., dF (x∗ )> Ψ = 0).
On s’intéresse à décrire l’allure de l’ensemble image F (U ) où U est un voisinage ouvert du point
x∗ . Bien sûr, si dF (x∗ ) est surjective, F (U ) est un voisinage ouvert de x∗ : c’est le théorème des
fonctions implicites.
Mais lorsque F (x∗ ) est au bord de F (IRn ) ce n’est pas le cas, et alors on va voir que F (U )
ressemble à l’intérieur d’une parabole, cf dessin.
72
F (U )
Ψ̃
F (x∗)
Im dF (x∗)
En effet, Im dF (x∗ ) décrit l’ensemble des points images par F , “à l’ordre 1" : c’est, de manière
approchée, l’ensemble des points F (x∗ + d) = F (x∗ ) + dF (x∗ ).d + o(kdk) où on néglige tous les
termes d’ordre > 1. Mais comme dF (x∗ ) n’est pas surjective, cet ensemble ne donne pas une bonne
approximation de l’ensemble image F (U ) ! (cf théorème des fonctions implicites) On a alors besoin
de faire une approximation d’ordre 2. Pour cela, on regarde la dérivée seconde intrinsèque, qui est
la forme quadratique :
Q = Ψ̃. d2 F (x∗ )| ker dF (x∗ )
C’est la Hessienne de F en x∗ , qu’on restreint au noyau de dF (x∗ ) (puisqu’on connaît déjà ce qui
se passe à l’ordre 1) et qu’on co-restreint, autrement dit qu’on regarde, le long de Ψ̃ (car, on le voit
sur le dessin, c’est la direction manquante). Du point de vue géométrie différentielle, cet objet est
intrinsèque : la forme quadratique Q est bien définie et ne dépend pas des choix de coordonnées.
Maintenant, on a les faits suivants :
— Si F (x∗ ) est au bord de F (IRn ) alors Q > 0 : la forme quadratique Q est positive.
— Si Q > 0, i.e., si la forme quadratique Q est définie positive, alors F (x∗ ) est “strictement au
bord" de F (U ), au sens où F (U ) ressemble à l’intérieur d’une parabole dont le sommet est
F (x∗ ) (cf dessin).
C’est une analyse d’ordre deux, de type “théorie de Morse" en géométrie différentielle, mais qui ne
fait que généraliser ce qu’on connaît bien sur les fonctions de IR dans IR : si x minimise f alors
f 0 (x) = 0 et f 00 (x) > 0, et réciproquement, si f 0 (x) = 0 et f 00 (x) > 0 alors x minimise f localement.
Il reste alors à relier Q à la Hessienne de L restreinte à ker dH(x∗ ). Tout d’abord, comme
F = (f, H) et que les coordonnées de Ψ̃ sont (1, λ∗1 , . . . , λ∗p , (µ∗j )j∈I(x∗ ) ), on a Ψ̃.F = L. Il reste
alors seulement à montrer que ker dF (x∗ ) = ker dH(x∗ ). Comme F = (f, H), on a ker dF (x∗ ) =
ker df (x∗ ) ∩ ker dH(x∗ ). Par ailleurs, x∗ vérifie la condition de multiplicateur de Lagrange, qui
s’écrit df (x∗ ) + Ψ̃.dH(x∗ ) = 0, ce qui implique que ker dH(x∗ ) ⊂ ker df (x∗ ), et donc ker dF (x∗ ) =
ker dH(x∗ ).
73
Le cadre est le suivant. On note s ∈ S le paramètre, où S est un espace de Banach. Pour tout
s ∈ S, on considère le problème d’optimisation sous contrainte d’égalité (on s’y ramène toujours
comme on l’a vu)
min {f (x, s) | x ∈ IRn t.q. h(x, s) = 0} (4.7)
où f : IRn × S → IR et h : IRn × S → IRp sont C 2 . Soit U ⊂ IRp un voisinage ouvert de 0. On
suppose ce problème “bien posé" au sens où, pour tout s ∈ S, il existe un unique minimiseur x(s).
On suppose aussi que, pour tout s ∈ S, le problème est qualifié, ce qui signifie que l’application
n p
linéaire ∂h
∂x (x(s), s) : IR → IR est surjective.
On définit la “fonction valeur"
On s’intéresse à savoir comment la fonction valeur S(s) et le minimiseur x(s) dépendent du para-
mètre s ∈ S : c’est ce qui s’appelle l’analyse de sensibilité (“sensitivity analysis" en anglais).
L’analyse de sensibilité se fait grâce au théorème des fonctions implicites appliqué au système
d’optimalité obtenu avec la règle des multiplicateurs de Lagrange. Pour tout s ∈ S, comme x(s)
est minimiseur, et comme le problème est supposé qualifié, il existe un multiplicateur de Lagrange
λ(s) ∈ IRp (vecteur colonne ; donc λ(s)> est un vecteur ligne) tel que
∂f ∂h
(x(s), s) + λ(s)> (x(s), s) = 0
∂x ∂x
En écriture Lagrangienne, on définit ici
p
X
L(x, λ, s) = f (x, s) + hλ, h(x, s)i = f (x, s) + λ> h(x, s) = f (x, s) + λi hi (x, s)
i=1
On voit maintenant qu’on est dans le cadre du théorème des fonctions implicites : on veut en effet
résoudre, par rapport à (x, λ), le système de n + p équations F (x, λ, s) = 0 à n + p inconnues (x, λ),
en fonction du paramètre s.
Pour appliquer le théorème des fonctions implicites, il faut vérifier que, pour tout s ∈ S, la
∂F
jacobienne ∂(x,λ) (x, λ, s) est inversible. Calculons cette jacobienne :
! ! !
∂2L ∂2L ∂2L ∂h > A C>
∂F
(x, λ, s) = ∂x2 (x, λ, s) ∂x∂λ (x, λ, s) = ∂x2 (x, λ, s) ∂x (x, s) =
∂(x, λ) ∂h ∂h ∂h C 0
∂x (x, s) ∂λ (x, s) ∂x (x, s) 0
74
2
avec A = ∂∂xL2 (x, λ, s) et C = ∂h
∂x (x, s).
Montrons que, sous la condition suffisante du second ordre du théorème 28, qui s’écrit ici
∂2L
A1 = (x, λ, s)| ker ∂h (x,s) définie positive
∂x2 ∂x
(plus généralement, il suffit de supposer que cette forme quadratique soit non dégénérée) la jaco-
∂F
bienne ∂(x,λ) (x, λ, s) est inversible. La preuve est similaire à celle qui a été faite en section [Link].
Soit (y, µ) ∈ IRn × IRp tel que ∂(x,λ)
∂F
(x, λ, s).(y, µ) = 0, i.e.,
Ay + C > µ = 0, Cy = 0.
Comme y ∈ ker C, la première équation s’écrit donc A1 y + C > µ = 0. On la multiplie par y > , et
comme y > C > = 0 on obtient y > A1 y = 0, d’où y = 0 puisque A1 est définie positive (non dégénérée
suffirait). Donc C > µ = 0, et comme C est supposée surjective, C > est injective et donc µ = 0.
On déduit du théorème des fonctions implicites que l’équation F (x, λ, s) = 0 se résout par
rapport à (x, λ) et donne s 7→ (x(s), λ(s)) de classe C 1 .
On a donc obtenu le résultat suivant.
Proposition 5. Soit s∗ ∈ S un paramètre fixé. On suppose que, pour s = s∗ , le problème (4.7)
admet un unique minimiseur local x∗ = x(s∗ ), ayant un multiplicateur de Lagrange (normal) λ(s∗ )
et que
∂h ∗ ∗ n p
— l’application linéaire ∂x (x(s ), s ) : IR → IR est surjective ;
∂2L ∗ ∗ ∗ ∗
— la forme quadratique ∂x2 (x (s ), λ(s ), s )| ker ∂h ∗ ∗
∂x (x(s ),s )
est définie positive.
∗
Alors il existe un voisinage ouvert V de s dans S tel que, pour tout s ∈ V , le problème (4.7)
admet un unique minimiseur local x(s), ayant un multiplicateur de Lagrange (normal) λ(s), tous
deux dépendant de s de manière C 1 . La fonction valeur S(s) dépend aussi de s de manière C 1 .
En dérivant la relation F (x(s), λ(s), s) = 0 par rapport à s ∈ S, on peut de plus noter que les
différentielles x0 (s) et λ0 (s) par rapport à s vérifient le système
! ! !
∂2L ∂h > x0 (s) ∂2L
∂x2 (x(s), λ(s), s) ∂x (x(s), s) =− ∂s∂x (x(s), λ(s), s)
∂h
∂x (x(s), s) 0 λ0 (s) ∂h
∂s (x(s), s)
ce qui est sous-jacent à la méthode numérique de résolution de Lagrange-Newton (voir plus loin).
Par ailleurs, en dérivant par rapport à s les relations S(s) = f (x(s), s) et h(x(s), s) = 0, on
obtient (on note toujours 0 la différentielle par rapport à s ∈ S)
∂f ∂f ∂h ∂h
S 0 (s) = (x(s), s).x0 (s) + (x(s), s) et (x(s), s).x0 (s) + (x(s), s) = 0.
∂x ∂s ∂x ∂s
Or, la relation de multiplicateur de Lagrange s’écrit ∂f > ∂h
∂x (x(s), s) + λ(s) ∂x (x(s), s) = 0, et en
0
l’appliquant à x (s) et en utilisant les relations ci-dessus, on obtient finalement
∂f ∂h
S 0 (s) = (x(s), s) + λ(s)> (x(s), s).
∂s ∂s
Remarque 46 (Cas particulier). Un cas particulier intéressant est lorsque f (x, s) = f (x) et
h(x, s) = h(x) − s avec S = IRp , autrement dit on considère le problème d’optimisation dont la
fonction valeur est
S(s) = min f (x).
h(x)=s
75
On peut le voir comme le problème d’optimisation minh=0 f où l’on perturbe les contraintes
d’égalité avec un petit paramètre s. Dans les conditions de la proposition 5, on obtient alors
S 0 (s) = −λ(s)> , c’est-à-dire
∂S
λi (s) = − (s) ∀i ∈ {1, . . . , p}.
∂si
Ce résultat est intéressant car il donne une interprétation des multiplicateurs de Lagrange : pour
i = 1, . . . , p, le multiplicateur λi = λi (0) associé au problème d’optimisation minh=0 f est égal à
∂S
λi = − ∂s i
(0), autrement dit il est lié à la façon dont la fonction valeur S est modifiée si on perturbe
la contrainte d’égalité hi (x) = 0 par hi (x) = si . En ce sens, le multiplicateur de Lagrange λi est
une mesure de la sensibilité de la fonction valeur par rapport à la contrainte hi (x) = 0.
Petit aparté : théorème de Danskin. Mentionnons ici un théorème parfois bien utile, qui
relève de l’analyse de sensibilité, qui permet de dériver une “fonction valeur". Le contexte est le
suivant. Soient X un espace de Banach, K 6= ∅ un espace topologique compact et f : X × K → IR
une fonction continue. On définit la fonction valeur S : X → IR par
Pour tout x ∈ X, on note K̂(x) = {y ∈ K | f (x, y) = S(x)} l’ensemble des minimiseurs de f (x, ·).
Il n’est pas difficile de montrer que la fonction S est continue.
Théorème 29 (Danskin). On suppose que f est Fréchet différentiable par rapport à sa première
variable et telle que sa différentielle de Fréchet ∂f
∂x par rapport à x soit continue sur X × K. Alors
la fonction S est différentiable dans toute direction au sens de Dini en tout point x ∈ X, et
0 S(x + th) − S(x) ∂f
S (x; h) = lim+ = min (x, y).h | y ∈ K̂(x) ∀h ∈ X.
t→0 t ∂x
De plus, étant donné x ∈ X, si K̂(x) = {yx } est un singleton (i.e., si f (x, ·) a un unique minimiseur
yx ∈ K), alors S est Fréchet différentiable en x et dS(x) = ∂f∂x (x, yx ).
Remarque 47. Si X = IR, K = [−1, 1] et f (x, y) = xy, alors S(x) = −|x|. En x = 0, S est
dérivable le long de toute direction au sens de Dini (elle est dérivable à droite et à gauche), mais
elle n’est pas directionnellement dérivable.
76
4.3 Algorithmes d’optimisation avec contraintes
Dans cette section nous allons étudier différentes méthodes algorithmiques pour les problèmes
d’optimisation avec contraintes. On distingue les méthodes primales (qui sont, en gros, de type
gradient) des méthodes duales (qui sont basées sur l’interprétation Lagrangienne).
min f (x)
x∈C
xk+1 = xk − ρk ∇f (xk ).
Le défaut est que, même si x0 ∈ C, les itérés suivants xk ne restent pas forcément dans C. L’idée de
la méthode de gradient projeté est toute simple : à chaque itération, on projette le nouveau point
dans le convexe C.
Soit PC la projection orthogonale sur C (voir le théorème du convexe dans le chapitre sur la
convexité : on rappelle que cette projection est bien définie et que, de plus, c’est une application
1-Lipschitzienne).
La méthode du gradient projeté est alors
et on obtient, comme on l’a vu précédemment, diverses variantes en prenant un pas fixe ou variable
(pas optimal par exemple).
Comme pour la méthode de gradient, on peut démontrer un théorème de convergence :
Théorème 30. On suppose que f est α-convexe (avec α > 0) et que ∇f est M -Lipschitzienne (avec
M > 0). 2 Alors le problème d’optimisation a une unique solution x∗ . Soient 0 < β1 < β2 < M 2α
2.
Si à toute itération on choisit le pas ρk ∈ [β1 , β2 ] alors, quel que soit le point initial x0 , la méthode
de gradient projeté converge vers x∗ .
Remarque 48. Comme pour la méthode de gradient, si les hypothèses sur f ne sont vraies que
localement autour de x∗ , alors la convergence est locale.
Démonstration. Tout d’abord, comme f est α-convexe, elle admet sur le convexe C un unique
minimiseur x∗ ∈ C. Par la condition nécessaire d’optimalité de f sur un convexe, on a
donc
hx∗ − t∇f (x∗ ) − x∗ , x − x∗ i 6 0 ∀x ∈ C ∀t > 0.
2. On peut noter que, forcément, α 6 M .
77
Par le théorème du convexe (caractérisation du projeté), on en déduit que
x∗ = PC (x∗ − t∇f (x∗ )) ∀t > 0. (4.8)
On a alors
xk+1 − x∗ = PC (xk − ρk ∇f (xk )) − PC (x∗ − ρk ∇f (x∗ ))
et comme PC est 1-Lipschitzienne, on en déduit que
kxk+1 − x∗ k2 6 kxk − x∗ − ρk (∇f (xk ) − ∇f (x∗ ))k2
6 kxk − x∗ k2 − 2ρk hxk − x∗ , ∇f (xk ) − ∇f (x∗ )i + ρ2k k∇f (xk ) − ∇f (x∗ )k2
Comme ∇f est M -Lipschitzienne, on a k∇f (xk ) − ∇f (x∗ )k 6 M kxk − x∗ k, et comme f est
α-convexe, on a hxk − x∗ , ∇f (xk ) − ∇f (x∗ )i > αkxk − x∗ k2 . Par conséquent,
kxk+1 − x∗ k2 6 (1 − 2αρk + M 2 ρ2k )kxk − x∗ k2 .
La fonction
2α
ϕ(t) = M 2 t2 − 2αt = M 2 t t − 2 ϕ(t)
M
α
2α β1 M2 β2 t
est négative et convexe sur l’intervalle [0, M 2 ], et s’annule
2α 2α
aux deux bords. Comme 0 < β1 6 ρk 6 β2 < M 2 , on en
M
déduit que
3. Voici une autre manière de calculer Pker(C) , le projecteur orthogonal sur ker(C) (avec C surjective) :
− Par définition du projecteur orthogonal, pour tout x, on a Pker(C) (x) ∈ ker C et x − Pker(C) (x) ⊥ ker C, i.e.,
CPker(C) (x) = 0 et ∀y, hx − Pker(C) (x), yi = 0.
− On a ker C = (Im C > )⊥ . En effet, x ∈ ker C ⇔ Cx = 0 ⇔ ∀y, hy, Cxi = 0 = hC > y, xi = 0 ⇔ x ∈ (Im C > )⊥ .
− Donc x − Pker(C) (x) ∈ (ker C)⊥ = Im C > , i.e., il existe y tel que x − Pker(C) (x) = C > y. Or, CPker(C) (x) = 0
donc Cx = CC > y d’où y = (CC > )−1 Cx (CC > est inversible car C est surjective), et donc Pker(C) (x) =
x − C > y = x − C > (CC > )−1 Cx.
78
• C = IRn+ : on a alors
où
sat(a, s, b) = min(max(a, s), b) s
a b
a si s 6 a 0
= s si a 6 s 6 b
b si b 6 s
(fonction de saturation).
Les contraintes de bornes sont très fréquentes en optimisation.
Remarque 50. Les contraintes de bornes ou de positivité sont d’autant plus importantes que, en
fait, tout problème d’optimisation sous contraintes d’égalité et d’inégalité non linéaires se ramène
à un problème d’optimisation en dimension plus grande, avec contraintes d’égalité et contraintes
de négativité sur les coordonnées. En effet, considérons le problème
avec h : IRn → IRp et g : IRn → IRq . Pour tout j ∈ {1, . . . , q}, la contrainte gj (x) 6 0 est
équivalente à gj (x) = yj et yj 6 0. Ainsi, en posant y = (y1 , . . . , yq ), le problème d’optimisation
(4.9) est équivalent à
min{f (x) | h(x) = 0, g(x) − y = 0, y 6 0}
qui est un problème d’optimisation en dimension n + q puisque, maintenant, l’inconnue est (x, y) ∈
IRn+q , comportant des contraintes d’égalité non linéaires, et des contraintes d’inégalité très simples :
négativité de certaines coordonnées.
Cette technique d’ajouter des “variables molles" (slack variables en anglais) est très utilisée en
optimisation.
Remarque 51. Au passage, décrivons une autre réduction, utilisant également des variables
molles. Le problème (4.9) se ramène à un problème avec contraintes d’égalité non linéaires seule-
ment. Pour cela, il suffit de dire que, pour tout j ∈ {1, . . . , q}, la contrainte gj (x) 6 0 est équivalente
à gj (x) = −yj2 , avec yj ∈ IR. Ainsi, en posant Gj (x, y) = gj (x)+yj2 et G = (G1 , . . . , Gq ), le problème
d’optimisation (4.9) est équivalent à
min f (x)
h(x)=0
G(x,y)=0
qui est un problème d’optimisation comportant uniquement des contraintes d’égalité non linéaires.
Cette technique peut d’ailleurs fournir une preuve alternative de KKT (mais, pour obtenir
le signe des multiplicateurs de Lagrange, il faut toutefois supposer toutes les fonctions de classe
C 2 , puis appliquer la condition nécessaire d’optimalité d’ordre 2 ; alors que le théorème KKT ne
nécessite que des fonctions de classe C 1 ).
79
Notons que, numériquement, cette technique, avec des yj2 qui sont non linéaires, n’est pas
forcément préférable à la technique précédente, qui donnait des contraintes de négativité très
simples.
Comme précédemment, le défaut est que, même si x0 ∈ C, les itérés suivants xk ne restent pas
forcément dans C. Comme précédemment, on projette chaque nouvel itéré sur C. On obtient la
méthode de Newton projetée :
Naturellement, on peut procéder de même avec les diverses variantes de la méthode de Newton
que nous avons vues : quasi-Newton, Barzilai-Borwein.
La remarque 49 s’applique aussi à ces différents cas.
où f : IR → IR est une fonction et C est un sous-ensemble de IRn (on ne suppose rien de particulier
n
ici). Notons χC : IRn → [0, +∞] la fonction (parfois appelée fonction indicatrice de C) définie par
0 si x ∈ C
χC (x) =
+∞ si x ∈ IRn \ C
Trivialement, le problème d’optimisation sous contraintes (4.10) est équivalent au problème d’op-
timisation sans contraintes
minn (f (x) + χC (x)) (4.11)
x∈IR
En effet, f (x) + χC (x) = +∞ dès que x ∈ / C, donc forcément le minimum est à chercher dans
l’ensemble C ; mais dès que x ∈ C, on a f (x) + χC (x) = f (x). D’où l’équivalence des deux problèmes
d’optimisation.
On a ainsi montré qu’un problème général d’optimisation sous contraintes est équivalent à un
problème d’optimisation sans contraintes ! Cependant, cela est au prix de manipuler la fonction
généralisée f + χC , qui prend ses valeurs dans IR ∪ {+∞}. En pratique cela ne résout donc rien car,
en général, le minimum est atteint au bord de C (on sature la contrainte ! si on ne la sature pas,
c’est alors qu’on avait un problème d’optimisation classique sans contraintes...), et c’est justement
au bord de C que la fonction f + χC n’est pas régulière : on ne peut la dériver en de tels points.
En pratique, l’idée est d’utiliser des fonctions ϕε : IRn → IR, indicées par ε > 0, qui “approchent"
la fonction indicatrice χC , i.e., pour tout x ∈ IRn , ϕε (x) → χC (x) lorsque ε → 0, et on considère la
famille (indicée par ε > 0) de problèmes pénalisés
80
On parle de pénalisation exacte dans le cas (4.11), lorsque la fonction de pénalisation est la
fonction indicatrice χC . On parle de pénalisation inexacte dans le cas (4.12), lorsqu’on utilise une
fonction de pénalisation ϕε qui approche la fonction indicatrice χC .
En supposant que le problème (4.10) a une unique solution x∗ et que, pour tout ε, le problème
pénalisé (4.12) a une unique solution xε , on s’attend à ce que xε → x∗ lorsque ε → 0. C’est bien
le cas sous des conditions naturelles. On a le résultat suivant.
Théorème 31. On suppose que f est continue et infinie à l’infini, et que C est fermé non vide.
Ainsi, le problème d’optimisation sous contraintes (4.10) a au moins une solution.
On suppose que, pour tout ε > 0, la fonction ϕε : IRn → [0, +∞) est continue et à valeurs
positives et que la famille de fonctions (ϕε )ε>0 converge vers la fonction indicatrice χC de C au
sens suivant :
• convergence simple sur IRn : pour tout x ∈ IRn , ϕε (x) → χC (x) quand ε → 0 ;
• convergence uniforme sur tout compact de l’ouvert IRn \ C :
Alors, pour tout ε > 0, le problème pénalisé (4.12) a au moins une solution xε ∈ IRn , la famille
(xε )ε>0 est bornée, et toute valeur d’adhérence de (xε )ε>0 est une solution optimale de (4.10). Si
de plus le problème (4.10) a une unique solution x∗ ∈ C alors xε → x∗ lorsque ε → 0.
Remarque 52. On a utilisé ici un seul paramètre de pénalisation ε, mais on aurait pu considérer
une famille à plusieurs paramètres. L’adaptation est triviale.
Démonstration. Comme f est continue et infinie à l’infini et C est fermé non vide, le problème
d’optimisation sous contraintes (4.10) a au moins une solution x∗ ∈ C (pas forcément unique).
Notons que χC (x∗ ) = 0 puisque x∗ ∈ C.
Montrons que la famille (xε )ε>0 est bornée dans IRn . Pour tout ε > 0, comme ϕε > 0, la
fonction f + ϕε est continue et infinie à l’infini, donc elle a au moins un minimiseur xε ∈ IRn (pas
forcément unique). En particulier, on a
où ε0 > 0 est fixé assez petit, car ϕε (x∗ ) → χC (x∗ ) = 0 quand ε → 0. Donc la famille de réels
(f (xε ))ε>0 est bornée, et comme f est continue et infinie à l’infini, on en déduit que la famille
(xε )ε>0 est bornée dans IRn .
Soit alors (xεk )k∈IN une sous-suite convergeant vers un point x̄ ∈ IRn . Montrons que x̄ est un
minimiseur du problème (4.10). D’après (4.13), on a f (xε ) 6 f (x∗ ) + ϕε (x∗ ), donc en passant à la
limite on obtient f (x̄) 6 f (x∗ ) = minf . Pour déduire que x̄ est un minimiseur du problème (4.10),
C
il manque le fait que x̄ ∈ C.
Montrons par l’absurde que x̄ ∈ C. Si x̄ ∈ / C : soit K une boule compacte de centre x̄, contenue
dans l’ouvert IRn \ C, et soit k0 ∈ IN assez grand pour que xεk ∈ K pour tout k > k0 . D’après
(4.13), on a f (xεk ) + ϕεk (xεk ) 6 f (x∗ ) + 1, et comme f (xεk ) → f (x̄) lorsque k → +∞, on en
déduit qu’il existe k1 > k0 tel que ϕεk (xεk ) 6 f (x∗ ) − f (x̄) + 2 pour tout k > k1 . Mais comme
xεk → x̄ lorsque k → +∞ et comme ϕεk converge uniformément vers +∞ sur le compact K, on a
ϕεk (xεk ) → +∞, d’où une contradiction avec l’inégalité précédente.
Par conséquent x̄ ∈ C. Comme f (x̄) 6 f (x∗ ), on doit avoir f (x̄) = f (x∗ ) (puisque x∗ minimise
f sur C) donc x̄ est solution du problème (4.10) (mais on n’a pas forcément x̄ = x∗ ).
Le théorème ci-dessus donne des conditions générales sur la fonction de pénalisation.
81
Une fois qu’on a choisi une fonction de pénalisation, on a un problème d’optimisation sans
contraintes, que l’on peut résoudre, par exemple, par une méthode de gradient (simple, pas optimal,
gradient conjugué) ou de Newton.
En pratique on utilise majoritairement les deux procédés suivants : pénalisation externe, ou
interne. On considère le problème d’optimisation
L’idée est ici que, lorsque ε > 0 est petit, dès que hi (x) 6= 0 ou que gj (x) > 0 alors ϕε (x) est très
grand. Donc, minimiser f + ϕε force hi (x) ' 0 et gj (x) . 0. On est bien dans le cadre du théorème
31. On parle ici de pénalisation externe, parce qu’un minimiseur xε de ϕε ne vérifie pas forcément
les contraintes. Il ne les vérifie que de manière approchée lorsque ε est petit.
Bien sûr, on peut imaginer de nombreuses variantes à la fonction ϕε ci-dessus. Par exemple
au lieu de choisir un seul paramètre de pénalisation ε > 0, on peut prendre différents paramètres
αi > 0, βj > 0, selon les indices des contraintes. On peut choisir aussi d’autres fonctions que les
carrés. Par exemple, on peut prendre
p q
X ai (hi (x)) X bj (gj (x))
ϕα,β (x) = +
i=1
αi j=1
βj
d
En effet, la fonction s 7→ max(s, 0)2 est dérivable sur IR, et ds max(s, 0)2 = 2 max(s, 0).
Pénalisation interne. Dans la pénalisation interne, il faut supposer que le point initial x0 de la
méthode itérative (par exemple gradient) qu’on met en oeuvre pour résoudre le problème pénalisé,
vérifie strictement les contraintes d’inégalité, i.e., est tel que g(x0 ) < 0.
Un exemple de fonction de pénalisation interne est
p q
1X 2
X 1
ϕε (x) = hi (x) + ε 2
. (4.17)
ε i=1 g
i=1 j
(x)
L’idée est que, au cours des itérations, si le point xk s’approche d’une frontière gj (x) = 0, alors
1
gj (xk )2 devient très grand, ce qui n’est pas favorisé par la minimisation de ϕε . On est de nouveau
82
dans le cadre du théorème 31. On parle ici de pénalisation interne, parce que les itérés xk , ainsi
que les minimiseurs de ϕε , restent à l’intérieur du domaine des contraintes en inégalité. Notons que
p q
X 2hi (x) X 2ε
∇ϕε (x) = ∇hi (x) − ∇gj (x). (4.18)
i=1
ε j=1
gj (x)3
Comme précédemment, de nombreuses variantes existent. Une variante très souvent utilisée est
la pénalisation logarithmique :
p q
1X X
ϕε (x) = hi (x)2 − ε ln(−gj (x)). (4.19)
ε i=1 j=1
C’est la même idée : si le point xk s’approche d’une frontière gj (x) = 0, alors − ln(−gj (x)) devient
très grand. Notons que
p q
X 2hi (x) X ε
∇ϕε (x) = ∇hi (x) − ∇gj (x). (4.20)
i=1
ε j=1
gj (x)
Remarque 53. Que ce soit pour une pénalisation externe ou interne, lorsqu’on met en oeuvre
une telle méthode de pénalisation, on doit :
1. choisir un paramètre de pénalisation ε > 0 ;
2. choisir une méthode de minimisation sans contrainte (par exemple, gradient simple ou pas op-
timal, gradient conjugué, Newton ou quasi-Newton), pour résoudre le problème de minimiser
f + ϕε sur IRn . Par exemple si on choisit une méthode de gradient simple, alors
Or, le choix de ε n’est en fait pas évident ! D’après le théorème de convergence 31, on a intérêt à
prendre ε > 0 petit. Mais si on prend ε trop petit, le calcul de ϕε (x) et de son gradient peuvent
générer des erreurs importantes (lorsqu’on divise par ε). En pratique, on utilise souvent une suite
(εk )k∈IN de paramètres de pénalisation, qui évolue au cours des itérations. Une méthode est de
diviser ε par 2, toutes les 10 itérations par exemple ; et si on observe que, à une itération donnée,
le problème est mal conditionné, alors on double la valeur de ε.
83
En comparant ces deux expressions, et en ayant en tête le théorème de convergence 31, on s’attend
donc à ce que
2hi (xε ) 2 max(gj (xε ), 0)
−→ λ∗i et −→ µ∗j
ε ε→0 ε ε→0
pour tous i, j. C’est le cas en effet, pourvu que le multiplicateur (normal) (λ∗ , µ∗ ) soit unique ! Cela
résulte d’un raisonnement facile par passage à la limite, calqué sur la deuxième démonstration de
KKT donnée en Section [Link].
Ce résultat est intéressant car il montre que la pénalisation donne un moyen d’estimer (de
manière approchée) les multiplicateurs de Lagrange.
Bien sûr, on a des résultats équivalents en utilisant d’autres pénalisations (cf les expressions
des gradients (4.18) ou (4.20)).
Remarque 54. En pratique, on utilise volontiers une méthode de pénalisation avec un ε pas
trop petit, pour avoir une première approximation (même grossière) du minimiseur x∗ et des
multiplicateurs de Lagrange (estimés comme ci-dessus). Cette estimation de (x∗ , λ, µ) peut alors
servir à initialiser, par exemple, une méthode de Lagrange-Newton (voir section suivante), plus
précise mais plus difficile à initialiser.
(dont la solution cherchée est (x∗ , λ∗ )) qui est un système de n + p équations à n + p inconnues
(x, λ).
La méthode de Lagrange-Newton consiste à résoudre le système d’équations (4.22) par la mé-
thode de Newton :
xk+1 xk
= − dF (xk , λk )−1 F (xk , λk )
λk+1 λk
Pour que ce problème de Newton soit bien posé (voir le théorème sur la méthode de Newton), il
faut assurer que la Jacobienne de F en (x∗ , λ∗ ) soit inversible, i.e., que la différentielle dF (x∗ , λ∗ ) :
IRn+p → IRn+p soit inversible. On sait alors que la méthode de Newton converge, dans un voisinage
de (x∗ , λ∗ ). Or, cette différentielle est
∂2L
∗ ∗ ∗ >
∗ ∗ 2 (x , λ ) dh(x )
dF (x , λ ) = ∂x
dh(x∗ ) 0
84
qui est la matrice de sensibilité vue en Section [Link] : elle est inversible sous les conditions que
2
dh(x∗ ) : IRn → IRp soit surjective et que la Hessienne ∂∂xL2 (x∗ , λ∗ ) restreinte à ker dh(x∗ ) soit définie
positive.
En appliquant le théorème de convergence de la méthode de Newton (théorème 22), on conclut
donc :
2
Proposition 6. Si dh(x∗ ) : IRn → IRp est surjective et si la Hessienne ∂∂xL2 (x∗ , λ∗ ) restreinte à
ker dh(x∗ ) est définie positive, alors il existe ε > 0 tel que, si k(x0 , λ0 )k 6 ε, alors l’algorithme de
Lagrange-Newton
∂2L −1
xk+1 xk (xk , λk ) dh(xk )> ∇x L(xk , λk )
= − ∂x2
λk+1 λk dh(xk ) 0 h(xk )
∂2L
sk = xk+1 − xk , gk = ∇f (xk ), hk = h(xk ), Ak = (xk , λk ), Ck = dh(xk ),
∂x2
le système linéaire (4.23) s’écrit
i.e.,
Ck>
Ak sk g
=− k (4.24)
Ck 0 λk+1 hk
| {z }
Mk
85
Variante : méthode de Wilson. Cette variante est basée sur la remarque suivante : le couple
(sk , λk+1 ) est solution du système (4.24) si et seulement si sk est la solution du problème d’opti-
misation sous contrainte d’égalité
1 > >
min s Ak s + gk s (4.25)
Ck s+hk =0 2
min L(x, λk )
h(x)=0
avec laquelle, à partir de l’itéré k, on déterminerait xk+1 . Ici, linéaire-quadratique veut dire :
approximation quadratique (i.e., à l’ordre 2) du critère de minimisation, et approximation linéaire
de la contrainte. En effet, connaissant xk , on cherche xk+1 = xk + sk tel que 0 ' h(xk+1 ) '
h(xk ) + dh(xk ).sk à l’ordre 1, ce qui donne la contrainte Ck sk + hk = 0. A l’ordre 2, on a
∂L 1 ∂2L
L(xk+1 , λk ) ' L(xk , λk ) + (xk , λk ).sk + (xk , λk ).(sk , sk )
∂x 2 ∂x2
et donc, dire que xk+1 minimise L par rapport à x, avec cette approximation d’ordre 2, revient à
dire que sk minimise
∂L 1 ∂2L
(xk , λk ).sk + (xk , λk ).(sk , sk ).
∂x 2 ∂x2
Or
∂2L
(xk , λk ).(sk , sk ) = s> k Ak sk ,
∂x2
et par ailleurs
∂L
(xk , λk ).sk = df (xk ).sk + λ> > > > >
k dh(xk ).sk = gk sk + λk Ck sk = gk sk − λk hk .
∂x
Ainsi, avec cette approximation, sk est bien le minimiseur de (4.25).
La méthode SQP (Sequential Quadratic Programming, en français : méthode de program-
mation quadratique séquentielle, ou successive) s’appuie sur cette idée d’approximation linéaire-
quadratique.
Considérons le problème d’optimisation sous contraintes d’égalité et d’inégalité
min f (x)
h(x)=0
g(x)60
86
Wilson ci-dessus, en faisant aussi une approximation linéaire des contraintes d’inégalité. L’itération
k s’écrit :
calculer sk , minimiseur de
2
1∂ L
min (xk , λk ).(s, s) + df (xk ).s
dh(xk ).s+h(xk )=0 2 ∂x2
dg(xk ).s+g(xk )60
avec les mêmes notations que d’habitude. Supposons que ce problème a un minimiseur x∗ , qui a
un multiplicateur de Lagrange normal λ∗ . Le Lagrangien de ce problème est alors
87
Toujours en notant que h(x∗ ) = 0, on calcule que la Hessienne par rapport à x du Lagrangien
augmenté au point (x∗ , λ∗ ) est
∂ 2 Lc ∗ ∗ ∂2L ∗ ∗
Qc = (x , λ ) = (x , λ ) + c dh(x∗ )> dh(x∗ ),
∂x2 ∂x2
i.e., c’est la somme dh(x∗ )> dh(x∗ ), qui est identifié à une forme quadratique symétrique positive,
2
et de ∂∂xL2 (x∗ , λ∗ ) qui, restreinte à ker dh(x∗ ), est supposée symétrique définie positive. Nous avons
le lemme suivant.
Lemme 5. Soient A ∈ Mn (IR) une matrice symétrique, et soit B ∈ Mp,n (IR), telles que si
x ∈ ker B \ {0} alors x> Ax > 0. Alors il existe c0 > 0 tel que pour tout c > c0 , A + cB > B est
symétrique définie positive.
Démonstration. Par l’absurde, supposons qu’il existe une suite (xk ) de IRn \ {0} telle que x>
k (A +
kB > B)xk 6 0. Quitte à diviser par kxk k2 , on peut supposer que kxk k = 1 pour tout k. Par
compacité, à sous suite-près, xk → x avec kxk = 1, et en passant à la limite dans l’inégalité
1 > > 2 > >
k xk (A + kB B)xk 6 0, on obtient kBxk = x B Bx 6 0 et donc Bx = 0, i.e., x ∈ ker B \ {0}.
Par ailleurs, on a xk Axk 6 −kxk B Bxk = −kkBxk k2 6 0 donc en passant à la limite, x> Ax 6 0.
> > >
On a donc obtenu l’existence d’un x ∈ ker B \{0} tel que x> Ax 6 0. Cela contredit l’hypothèse.
2
D’après ce lemme, il existe c0 > 0 tel que, pour tout c > c0 , la Hessienne ∂∂xL2c (x∗ , λ∗ ) est
symétrique définie positive. Comme ∇x Lc (x∗ , λ∗ ) = 0, il s’ensuit que x∗ est un minimiseur local
strict de la fonction Lc (·, λ∗ ).
Remarque 57. Dans les conditions du théorème, on a la propriété Lc (x∗ , λ∗ ) 6 Lc (x, λ∗ ) pour
tout x ∈ V où V est un voisinage de x∗ . De plus, comme h(x∗ ) = 0, on a aussi Lc (x∗ , λ∗ ) = Lc (x∗ , λ)
pour tout λ ∈ IRp . On a donc la double inégalité
qui signifie que (x∗ , λ∗ ) est un point selle (local) du Lagrangien augmenté Lc (voir section [Link]
plus loin). La propriété de point selle est à la base des méthodes duales étudiées en section 4.3.2.
Déduisons de cette propriété l’algorithme du Lagrangien augmenté. Supposons que, à l’itération
k, on dispose d’un paramètre de pénalisation ck et d’un multiplicateur λk . Le théorème (32) suggère
de minimiser la fonction x 7→ Lck (x, λk ), donc, on cherche
ce qui suggère de choisir λk+1 = λk + ck h(xk ), car en effet, on cherche à faire en sorte que xk ' x∗
et λk+1 ' λ∗ , vérifiant la condition de multiplicateurs de Lagrange df (x∗ ) + (λ∗ )> dh(x∗ ) = 0.
Finalement, l’algorithme du Lagrangien augmenté est :
Cet algorithme est assez proche de l’algorithme d’Uzawa dans les méthodes duales (voir la section
[Link]). Il existe de nombreuses variantes. Tout d’abord, dans l’algorithme ci-dessus, on peut
88
choisir la suite (ck ) de diverses manières, par exemple, constante assez grande ; ou bien croissante
et majorée. Ensuite, le problème de minimisation en x peut être résolu de beaucoup de manières.
Si on se contente d’une approximation, on peut par exemple faire un pas de gradient (de pas donné
ρ) par rapport à l’itération précédente :
où f : IRn → IR, g : IRn → IRq et h : IRn → IRp . Comme on supposera dans la suite qu’on est dans
le cas normal, on définit le Lagrangien
Remarque 58. On suppose qu’il existe un minimiseur x∗ du problème (4.28), associé à un mul-
tiplicateur normal (λ∗ , µ∗ ) ∈ IRp × IRq+ (i.e, λ0 = 1) vérifiant la condition de complémentarité
(µ∗ )> g(x∗ ) = 0. On a h(x∗ ) = 0, et comme g(x∗ ) 6 0, on a aussi µ> g(x∗ ) 6 0 pour tout µ ∈ IRq+ ,
et donc
L(x∗ , λ, µ) = f (x∗ ) + λ> h(x∗ ) + µ> g(x∗ ) 6 f (x∗ ) = L(x∗ , λ∗ , µ∗ ) ∀(λ, µ) ∈ IRp × IRq+
(pour la dernière égalité L(x∗ , λ∗ , µ∗ ) = f (x∗ ), on a utilisé les égalités h(x∗ ) = 0 et (µ∗ )> g(x∗ ) = 0).
Cette inégalité montre que (λ∗ , µ∗ ) maximise la fonction (λ, µ) 7→ L(x∗ , λ, µ) sur IRp × IRq+ . On
avait déjà vu ce fait dans la remarque 43.
Par ailleurs, dans cette même remarque on a dit que x∗ ne minimise pas forcément la fonction
x 7→ L(x, λ∗ , µ∗ ) sur IRn , sauf si f et les fonctions gj sont convexes et les fonctions hi sont affines
(voir proposition 4), ou bien, sauf si on utilise le Lagrangien augmenté (voir section [Link]) sous
condition de second ordre.
Cette remarque préliminaire motive la définition de point selle de L.
89
On appelle ce genre de point, un point selle (saddle point
en anglais) ou un point col.
A droite, on voit le graphe de la fonction
f (x, y) = x2 − y 2 .
Théorème 33. Si (x∗ , λ∗ , µ∗ ) ∈ IRn × IRp × IRq+ est un point selle de L alors x∗ est un minimiseur
du problème (4.28). Si de plus f , g et h sont différentiables, ce minimiseur x∗ vérifie les conditions
KKT avec le multiplicateur normal (λ∗ , µ∗ ).
Démonstration. Montrons d’abord que x∗ vérifie les contraintes, i.e., montrons que h(x∗ ) = 0 et
g(x∗ ) 6 0. Comme L(x∗ , λ, µ) 6 L(x∗ , λ∗ , µ∗ ) pour tout (λ, µ) ∈ IRp × IRq+ , on a
En particulier, pour x ∈ IRn vérifiant les contraintes h(x) = 0 et g(x) 6 0 (donc (µ∗ )> g(x) 6 0),
on en déduit que f (x∗ ) 6 f (x), ce qui est le résultat souhaité.
Lorsque f , g et h sont différentiables, le fait que x∗ minimise L(·, λ∗ , µ∗ ) sur IRn implique que
∇x L(x∗ , λ∗ , µ∗ ) = 0, ce qui donne KKT avec le multiplicateur normal (λ∗ , µ∗ ).
90
Ce résultat complète la proposition 4 et la remarque préliminaire 58.
Démonstration. Notons au préalable que, sous ces hypothèses, le Lagrangien L(x, λ, µ) = f (x) +
λ> h(x)+µ> g(x) est convexe par rapport à x (car µ ∈ IRq+ ). Montrons que (ii) ⇒ (iii) ⇒ (i) ⇒ (ii).
Partons de (ii). Comme x∗ est un minimiseur du problème (4.28) (qui est qualifié), alors il
existe (λ∗ , µ∗ ) ∈ IRp × IRq+ tel que ∇x L(x∗ , λ∗ , µ∗ ) = 0 (égalité des multiplicateurs de Lagrange)
et µ∗j gj (x∗ ) = 0 pour tout j ∈ {1, . . . , q} (conditions de complémentarité). Ainsi, (x∗ , λ∗ , µ∗ )
vérifie KKT, ce qui est (iii). Notons que l’égalité h(x∗ ) = 0 et les conditions de complémentarité
impliquent que L(x∗ , λ∗ , µ∗ ) = f (x∗ ). Par convexité de L, la condition ∇x L(x∗ , λ∗ , µ∗ ) = 0 est
équivalente à dire que x∗ minimise L(·, λ∗ , µ∗ ), i.e.,
Par ailleurs, d’après la remarque préliminaire 58, (λ∗ , µ∗ ) maximise la fonction (λ, µ) 7→ L(x∗ , λ, µ)
sur IRp ×IRq+ . Donc (x∗ , λ∗ , µ∗ ) est un point selle de L, ce qui est (i). Enfin, (i) ⇒ (ii) par le théorème
33.
Bien que les hypothèses de la proposition 7 ne soient pas vérifiées en général, comme on l’a
souvent fait précédemment, on considère qu’elles sont vraies au moins localement et approximati-
vement, ce qui nous incite donc à chercher des points selles de L, i.e., à minimiser L par rapport
à x ∈ IRn (minimisation sans contrainte) et maximiser L par rapport à (λ, µ) ∈ IRp × IRq+ (maxi-
misation sous contrainte de positivité).
C’est le principe des méthodes duales qu’on va voir plus loin.
D’autre part, bien que la réciproque du théorème 33 ne soit pas vraie en général, on a vu en
section [Link] que, en ramenant au préalable le problème comme un problème d’optimisation sous
contraintes d’égalité seulement, on peut considérer le problème équivalent (4.27) qui conduit à
définir le Lagrangien augmenté Lc , et on a vu en remarque 57 que, dans les conditions du théorème
32, le Lagrangien augmenté Lc admet un point selle (local) qui est exactement (x∗ , λ∗ ). On peut
donc en fait toujours se ramener à cette situation, en tout cas au moins localement.
et que le problème d’optimisation sous contraintes (4.28) est équivalent au problème de minimiser
f + χC sur IRn (problème sans contrainte), i.e.,
91
En effet, dès que h(x) 6= 0 on a supλ∈IRp λ> h(x) = +∞ (et sinon il vaut 0), et dès que gj (x) > 0
pour un j ∈ {1, . . . , q} on a supµ∈IRq+ µ> g(x) = +∞ (et sinon il vaut 0). On peut donc écrire le
problème (4.28) sous la forme
i.e., en reconnaissant le Lagrangien L(x, λ, µ) = f (x) + λ> h(x) + µ> g(x), et en écrivant un inf
plutôt qu’un min,
inf f (x) = inf n sup L(x, λ, µ) (4.29)
h(x)=0 x∈IR λ∈IRp
g(x)60 q
µ∈IR
+
Définition 10. Le problème (4.29) s’appelle le problème primal. Il coïncide avec le problème
d’optimisation initial (4.28).
On observe que le problème primal s’écrit comme “inf sup L" : on maximise d’abord L par
rapport à (λ, µ), puis on minimise par rapport à x.
On pourrait fort bien considérer de faire l’inverse ! Cela s’appelle alors le problème dual :
Définition 11. On appelle problème dual le problème
Dans la définition du problème dual, on minimise d’abord L par rapport à x, puis on maximise
par rapport à (λ, µ). On a toujours
autrement dit, la valeur optimale du problème dual est toujours inférieure ou égale à la valeur
optimale du problème primal. En effet, pour tout (x0 , λ, µ) ∈ IRn × IRp × IRq+ , on a
On prend alors l’inf sur les x0 ∈ IRn dans cette inégalité, et on obtient (4.31).
Saut de dualité. On n’a pas forcément égalité dans (4.31). La différence s’appelle saut de dualité.
On va voir ci-dessous un théorème assurant que le saut de dualité est nul dans certaines conditions.
Notons que, même lorsqu’il y a un saut de dualité non nul, résoudre le problème dual est tout
de même intéressant car sa valeur optimale fournit un minorant de la valeur optimale du problème
primal.
92
[Link] Théorème de dualité
La fonction
w(λ, µ) = inf n L(x, λ, µ)
x∈IR
s’appelle fonction duale. Le problème dual consiste à maximiser w(λ, µ) sur IRp × IRq+ . Pour tout
x ∈ IRn fixé, la fonction (λ, µ) 7→ L(x, λ, µ) est affine, donc w est un infimum de fonctions affines,
donc w est concave (on a déjà vu ce fait).
Par contre, w n’est pas forcément différentiable. Toutefois, pour un (λ, µ) ∈ IRp × IRq+ fixé,
si x est l’unique minimiseur de L(·, λ, µ) sur IRn , alors d’après le théorème 29 de Danskin, w est
différentiable en (λ, µ), et ∇λ w(λ, µ) = h(x) et ∇µ w(λ, µ) = g(x).
Théorème 34. S’il existe un point selle (x∗ , λ∗ , µ∗ ) ∈ IRn × IRp × IRq+ de L, alors le saut de dualité
est nul, x∗ est solution de (4.28) (associé au multiplicateur de Lagrange normal (λ∗ , µ∗ ) si de plus
f , g et h sont différentiables) et
sup inf n L(x, λ, µ) = inf n sup L(x, λ, µ) = L(x∗ , λ∗ , µ∗ ) = f (x∗ ) = w(λ∗ , µ∗ ) = sup w(λ, µ)
λ∈IRp x∈IR x∈IR λ∈IRp λ∈IRp
q q q
µ∈IR µ∈IR µ∈IR
+ + +
Réciproquement, si (4.28) a un minimiseur x∗ et s’il existe (λ∗ , µ∗ ) ∈ IRp × IRq+ tel que f (x∗ ) =
w(λ∗ , µ∗ ) alors (x∗ , λ∗ , µ∗ ) est un point selle de L.
Démonstration. Supposons que (x∗ , λ∗ , µ∗ ) ∈ IRn × IRp × IRq+ est un point selle de L. D’après le
théorème 33, on sait déjà que x∗ est un minimiseur de (4.28), et on a vu dans la preuve de ce
théorème que, de plus, µ∗j gj (x∗ ) = 0 pour tout j ∈ {1, . . . , q} (conditions de complémentarité),
donc en particulier (µ∗ )> g(x∗ ) = 0, et f (x∗ ) = L(x∗ , λ∗ , µ∗ ). Comme (x∗ , λ∗ , µ∗ ) est un point selle
de L, on a
f (x∗ ) = L(x∗ , λ∗ , µ∗ ) = inf n L(x, λ∗ , µ∗ ) = w(λ∗ , µ∗ ).
x∈IR
et donc on n’a que des égalités ci-dessus, ce qui est le résultat cherché.
Réciproquement, supposons que (4.28) a un minimiseur x∗ et qu’il existe (λ∗ , µ∗ ) ∈ IRp × IRq+
tel que f (x∗ ) = w(λ∗ , µ∗ ). Par définition,
(car h(x∗ ) = 0), donc 0 6 (µ∗ )> g(x∗ ), mais comme µ∗ > 0 et g(x∗ ) 6 0, on a forcément
µ∗j gj (x∗ ) = 0 pour tout j ∈ {1, . . . , q} (conditions de complémentarité). Donc L(x∗ , λ∗ , µ∗ ) = f (x∗ ).
De nouveau, par hypothèse et par définition, f (x∗ ) = w(λ∗ , µ∗ ) 6 L(x, λ∗ , µ∗ ) pour tout x ∈ IRn ,
donc L(x∗ , λ∗ , µ∗ ) 6 L(x, λ∗ , µ∗ ) pour tout x ∈ IRn , ce qui est l’une des propriétés de point selle.
D’autre part, comme h(x∗ ) = 0 et g(x∗ ) 6 0, on a
pour tout (λ, µ) ∈ IRp × IRq+ , donc L(x∗ , λ, µ) 6 L(x∗ , λ∗ , µ∗ ) pour tout (λ, µ) ∈ IRp × IRq+ . On a
démontré la propriété de point selle.
93
C’est la première partie du théorème 34 qui est la plus importante. Elle montre que, lorsqu’on a
un point selle de L, pour déterminer le triplet optimal (x∗ , λ∗ , µ∗ ) on peut résoudre indifféremment
le problème inf sup L ou sup inf L : on peut, dans n’importe quel ordre, minimiser L par rapport à
x ∈ IRn et maximiser L par rapport à (λ, µ) ∈ IRp × IRq+ .
La première étape peut avoir, naturellement, beaucoup de variantes. On peut utiliser n’importe
quelle méthode de minimisation sans contraintes. Toutefois, cette première étape nécessite de faire,
potentiellement, un grand nombre d’itérations. La variante suivante est alors très naturelle.
xk+1 = xk − ρk ∇x L(xk , λk , µk )
λk+1 = λk + ρk ∇λ L(xk , λk , µk ) = λk + ρk h(xk )
µk+1 = max(µk + ρk ∇µ L(xk , λk , µk ), 0) = max(µk + ρk g(xk ), 0)
94
Chapitre 5
Conclusion et compléments
95
Les méthodes d’optimisation globale sont souvent basées sur des heuristiques (on parle de
“métaheuristique"), inspirées par des systèmes naturels, par exemple :
— en physique : méthode du recuit simulé, s’inspirant de processus utilisés en métallurgie ;
— en biologie évolutionnaire : algorithmes génétiques ;
— en éthologie (étude du comportement des espèces) : algorithmes de colonies de fourmis,
en anglais “ant colony optimization", ou bien la méthode d’optimisation par essaims
particulaires, en anglais “particle swarm optimization".
Ces dernières méthodes connaissent actuellement un développement important, en lien avec
les études de dynamiques collectives (alignement, consensus, auto-organisation). Dans le do-
maine de l’Intelligence Artificielle (IA), on parle même de “swarm intelligence" (intelligence
de l’essaim), faisant référence à la capacité d’un groupe, d’un système connecté, a priori
décentralisé, à s’auto-organiser.
• On a étudié des problèmes d’optimisation avec un seul critère à minimiser (ou maximiser).
Dans la théorie de l’optimisation multi-objectifs, ou multi-critères, on a plusieurs fonctions à
minimiser. Comme on ne peut généralement pas trouver de point qui les minimise toutes en
même temps, on définit des ordres de préférence sur les objectifs, ce qui conduit à différentes
stratégies possibles : équilibres de Pareto, de Nash, de Stackelberg (liés à la théorie des jeux).
Naturellement, tous les points mentionnés ci-dessus peuvent être combinés.
Par ailleurs, on souhaite appliquer l’optimisation (notamment) à des problèmes en grande di-
mension : science des données, data science, big data, problèmes liés à l’IA. Dans de tels problèmes,
outre les difficultés et nouveautés mentionnées ci-dessus, une autre difficulté surgit rapidement :
celle de calculer les différentielles et/ou Hessiennes des fonctions en jeu. Jusqu’à présent, dans les
méthodes de gradient, de Newton, Uzawa, etc, on a supposé qu’on savait calculer explicitement les
différentielles et Hessiennes des fonctions. Mais en pratique, ce calcul peut s’avérer difficile, ou être
une source importante d’erreurs de codage (erreurs d’indices par exemple, lorsque les fonctions sont
données de manière discrétisée et qu’on veut faire des différences finies). Pour pallier à ce problème
on peut combiner des routines d’optimisation à des routines de différentiation automatique. C’est
l’objet de la section 5.1 ci-dessous.
f (x + h) − f (x)
f 0 (x) = + o(1)
h
96
lorsque h → 0, qui repose, comme on le sait, sur le développement limité à l’ordre 1
on suppose ici que f admet une extension holomorphe (au moins au voisinage de x) et on écrit le
développement limité à l’ordre 2, en complexes,
h2 00
f (x + ih) = f (x) + ihf 0 (x) − f (x) + o(h2 ).
2
La grosse différence maintenant est que
f (x + ih)
f 0 (x) = Im + o(h)
h
On voit qu’on a gagné un ordre dans l’approximation ! De plus, ce n’est plus une différence finie :
dans la formule ci-dessus, on prend la partie imaginaire de f (x + ih) (à laquelle, pour des raisons
informatiques, on peut accéder directement) qu’on divise par h ; alors que, dans la différence finie,
on faisait la différence f (x + h) − f (x) : autrement dit, la somme de deux termes, ce qui génère
d’importantes erreurs d’arrondi. Ainsi, la différentiation automatique, basée sur ce principe, permet
de calculer des dérivées à l’ordre de précision de la machine, soit à 10−14 près.
Cette idée (remarquable) a été largement développée depuis les années 60. Des dizaines d’années
de recherches ont about à des outils sophistiqués, très efficaces, pour dériver des fonctions (la
difficulté principale étant de calculer des dérivées de fonctions composées, mais ce n’est pas le lieu
de développer cet aspect ici). AMPL inclut de tels procédés de différentiation automatique. C’est
très intéressant car, dans un code, cela évite d’avoir à fournir les différentielles et Hessiennes des
fonctions en jeu.
Comme dit ci-dessus, AMPL permet de faire appel à une routine d’optimisation de notre choix, et
calcule par différentiation automatique les différentes dérivées nécessaires. Il existe de nombreuses
routines d’optimisation, très efficaces (et qui résultent de dizaines d’années de développements),
dans des domaines divers. En optimisation non linéaire sous contraintes non linéaires, à l’heure
actuelle, les deux routines les plus efficaces semblent être IpOpt et Knitro.
En gros, IpOpt (Interior Point Optimizer ) est une routine d’optimisation (désormais codée en
C++), basée sur une méthode de point intérieur (pénalisation interne ; mais c’est une routine très
élaborée), voir la documentation :
[Link]
Elle est gratuite et peut être installée sur n’importe quelle machine, voir
[Link]
Knitro (Nonlinear Interior point Trust Region Optimization) est une routine commerciale
(payante) qui, majoritairement, est basée sur une méthode SQP :
[Link]
Du point de vue efficacité générale sur des problèmes généraux, IpOpt et Knitro ont des
performances comparables. Actuellement, elles sont considérées comme étant les plus efficaces pour
résoudre des problèmes généraux d’optimisation non linéaire sous contraintes non linéaires.
Utilisation pratique. On peut utiliser AMPL combiné avec une routine d’optimisation (on conseille
IpOpt qui est gratuit) pour résoudre très facilement des problèmes d’optimisation non linéaire. Co-
der en AMPL est très simple, voici un exemple qui se passe de commentaire :
97
var x1; #, default 0;
var x2; #, default 0;
# Choix du solver :
option solver ipopt; # ligne a commenter si on lance le code sur le site AMPL
solve;
98
Apprendre à coder en AMPL n’est vraiment pas difficile. On trouve sur le web de nombreux
exemples de routines AMPL, par exemple :
[Link]
[Link]
par exemple [Link]
[Link]
et beaucoup d’autres. En particulier, on peut étudier le livre d’AMPL,
[Link]
99
printf: " %24.16e %24.16e\n", T, u[Nt-1] >> [Link];
printf{i in 0..Nt, j in 0..Nx}: " %24.16e\n", y[i,j] >> [Link];
end;
Ci-dessus, on a choisi de renvoyer les résultats dans un fichier texte [Link], qu’on lit ensuite
dans Matlab pour l’affichage graphique, à l’aide du fichier Matlab ci-dessous :
tempscont = tempscont’;
t=tempscont(:,1); u=tempscont(:,2);
y = reshape(mat_y,Nx+1,Nt+1); y=y’;
x = 0:1/Nx:1;
Bien entendu, on peut procéder différemment. L’affichage graphique des résultats peut se faire
dans Scilab, ou bien directement avec Python.
100
Beaucoup d’astuces et de manières de coder existent, on trouve quantité de matériel et de
“templates" sur le web, dont on pourra s’inspirer.
AMPL est un langage d’une très grande simplicité et d’une très grande puissance, mais sa licence
est payante. Il existe toutefois d’autres alternatives, gratuites, pour faire de la différentiation auto-
matique. Les solutions existantes sont nombreuses mais en général n’atteignent pas l’efficacité ou
la simplicité d’AMPL, à l’exception notable de CasADi :
[Link]
qui est une excellente solution de différentiation automatique pour faire de l’optimisation et du
contrôle optimal, en étant combinable à une routine d’optimisation comme IpOpt. Bien qu’étant
un peu moins simple d’utilisation que AMPL, CasADi est aussi efficace et est utilisable en Python
(ou même directement en C++). Le site web contient de nombreux exemples d’utilisation.
autrement dit, pour tout x ∈ IRn , F (x) est l’espérance de la variable aléatoire f (x, W ), (Ω, P ) est
un espace probabilisé (univers d’évènements ω), W : Ω → X est une variable aléatoire, et µ = W∗ P
est la loi de W (loi de probabilité sur X qui est l’image de P par W ).
On suppose f différentiable par rapport à x, de sorte que
Z Z
∇F (x) = E∇x f (x, ·) = ∇x f (x, W (ω)) dP (ω) = ∇x f (x, w) dµ(w).
Ω X
En pratique, pour estimer une espérance, on utilise une méthode de Monte-Carlo et on procède à
un échantillonnage, i.e., on fait une sélection de N évènements ωi ∈ Ω qui donnent des échantillons
wi = W (ωi ), i = 1, . . . , N (SAA : Sample Average Approximation). Lorsque µ est la mesure de
probabilité uniforme sur X, on a alors (en approximation)
N N
1 X 1 X
F (x) = fi (x) et ∇F (x) = ∇fi (x)
N i=1 N i=1
avec fi (x) = f (x, wi ). Autrement dit, la fonction F à minimiser est la moyenne d’un très grand
nombre N de fonctions fi . En machine learning, N est le nombre de données à partir desquelles
on veut “apprendre" les paramètres x ∈ IRn du modèle (voir section 5.3).
Pour minimiser la fonction F , un algorithme de descente de gradient à pas variable tk > 0
serait du type
xk+1 = xk − tk ∇F (xk )
mais le calcul de ∇F (xk ) est lourd, puisqu’il réclame de calculer la moyenne de tous les gradients
∇x f (x, w), c’est-à-dire, dans le cas échantillonné, le calcul des N gradients ∇fi (xk ), alors que N
est très grand.
101
L’algorithme du gradient stochastique n’utilise qu’un seul gradient pour chaque itération :
où wk+1 = W (ωk+1 ) ∈ X est un évènement aléatoire, tiré aléatoirement dans l’espace probabilisé
(X, µ), indépendamment du tirage précédent νk . Les tirages aléatoires wk = W (ωk ) forment une
suite (wk )k∈IN∗ qui est une réalisation d’un échantillon (de taille infinie) de la suite de variables
aléatoires (Wk )k∈IN∗ i.i.d. (“independently identically distributed"), c’est-à-dire des variables aléa-
toires indépendantes et qui ont toutes la même loi que W . Dans le cas échantillonné mentionné
ci-dessus, cela s’écrit
xk+1 = xk − tk ∇fik+1 (xk )
où ik+1 est un indice tiré aléatoirement et uniformément dans {1, . . . , N }, chaque tirage d’indice
étant indépendant du précédent.
Théorème 35. On suppose que F est C 1 et α-convexe : il existe donc un unique minimiseur
x∗ ∈ IRn . On suppose qu’il existe C > 0 telle que
Z
Ek∇x f (x, ·)k2 = k∇x f (x, w)k2 dµ(w) 6 C + Ckx − x∗ k2 ∀x ∈ IRn . (5.2)
X
(par exemple, tk = 1
k+1 ) alors la suite (xk )k∈IN converge presque sûrement vers x∗ .
Dans le cas échantillonné, l’hypothèse (5.2) s’écrit
N
1 X
k∇fi (x)k2 6 C + Ckx − x∗ k2 ∀x ∈ IRn
N i=1
avec C indépendant de N . Dans la littérature, il est souvent supposé qu’il existe C > 0 tel que
kxk+1 − x∗ k2 = kxk − x∗ k2 − 2tk hxk − x∗ , ∇x f (xk , wk+1 )i + t2k k∇x f (xk , wk+1 )k2 .
Prenons
R l’espérance de cette égalité par rapport à wk+1 : autrement dit, on applique l’intégration
X
(·) dµ(w k+1 ) à cette égalité. Par hypothèse d’indépendance, xk ne dépend pas de wk+1 , donc
Ekxk+1 − x∗ k2 = kxk − x∗ k2 − 2tk hxk − x∗ , E∇x f (xk , ·)i + t2k Ek∇x f (xk , ·)k2
| {z }
∇F (xk )
d’une part, comme F est α-convexe, on a h∇F (xk ) − ∇F (x∗ ), xk − x∗ i > αkxk − x∗ k2 avec
∇F (x∗ ) = 0 puisque x∗ minimise F . D’autre part, on utilise l’hypothèse (5.2), et on en déduit que
102
En prenant maintenant
R l’espérance de cette inégalité par rapport à wk (autrement dit, en appliquant
l’intégration X (·) dµ(wk ) à cette inégalité), et en posant
lorsque k → +∞, ce qui est vérifié sous les hypothèses sur la suite de pas. En effet, on a
k
X
ln Pk = ln(1 − 2αtj + Ct2j )
j=0
Comme tj → 0 lorsque j → +∞, on a ln(1 − 2αtj + Ct2j ) = −2αtj + O(t2j ) lorsque j → +∞, et en
sommant, en utilisant le fait que tj > 0, que la série des tj diverge et que la série des t2j converge,
Pk
on obtient ln Pk = −2α j=0 tj + O(1) lorsque k → +∞ et donc il existe des constantes C1 > 0 et
C2 > 0 telles que
Xk Xk
C1 exp − 2α tj 6 Pk 6 C2 exp − 2α tj ∀k ∈ IN
j=0 j=0
1
Pour simplifier, montrons la convergence vers 0 du membre de droite dans le cas où ti = i+1 (ce
n’est déjà pas si facile !). On sait que
k
X 1
= ln k + γ + o(1)
j=1
j
et donc Cst 1
k2α si α< 2
k k k
X X Cst X
t2i exp − 2α tj ∼ i2α−2 ∼ Cst lnkk si α= 1
2
i=0 j=i+1
k 2α i=1
Cst 1
si α>
k 2
103
car, par comparaison entre série et intégrale, on a
1
Cst si α< 2
X k
2α−2 1
i ∼ Cst ln k si α= 2
i=1 1
Cst k 2α−1 si α>
2
Dans tous les cas, le terme converge vers 0, ce qui est la conclusion désirée.
A ce stade, on a donc montré que (xk )k∈IN converge en moyenne quadratique vers x∗ , i.e.,
Ekxk − x∗ k2 → 0 lorsque k → +∞. Pour montrer que xk → x∗ presque sûrement, on utilise un
résultat général sur les martingales en théorie des probabilités. On pose
+∞
X
Zk = kxk − x∗ k2 + C t2i
i=k
donc (Zk )k∈IN est une sur-martingale. Comme elle est minorée (car positive), elle converge presque
sûrement, ce qui achève la démontration du théorème.
On rappelle ici qu’une martingale est une suite (Zk )k∈IN de variables aléatoires qui vérifie
E(Zk+1 | Z0 , . . . , Zk ) = Zk . On parle de sur-martingale si l’égalité est remplacée par 6, et de
sous-martingale si l’égalité est remplacée par >. Il est connu, en théorie des probabilités, que
toute sous-martingale majorée converge presque sûrement ; toute sur-martingale minorée converge
presque sûrement ; toute martingale majorée ou minorée converge presque sûrement.
On peut noter toutefois que la convergence de l’algorithme de gradient stochastique est lente.
Il faut donc trouver un compromis entre cette lenteur d’exécution et le nombre N de données.
Lorsque N est grand, le coût d’une itération de gradient stochastique est N fois plus petit que le
coût d’une itération de gradient classique (dans lequel on calculerait N gradients).
ou de manière plus générale, on se donne une fonction coût (souvent appelée “loss function" en
anglais) et on considère le problème
min C(f (X), Y )
f ∈F
f (X) = σk Ak · · · σ2 A2 σ1 A1 X
104
où Ai est une matrice de taille ni × ni−1 (autrement dit, Ai : IRni−1 → IRni est une application
linéaire), et σi : IRni → IRni est une application non linéaire, pour i = 1, . . . , k. Ci-dessus, pour
éviter la confusion de parenthèses, la notation σi a signifie σi (a).
L’entier k est appelé le nombre de couches (“layers" en anglais), les coefficients des matrices
Ai sont appelés poids (ils représentent généralement des interactions dans un réseau de neurones),
et les applications non linéaires σi sont appelées fonctions d’activation. Généralement, pour a =
(a1 , . . . , ani ) ∈ IRni , on prend σi (a) = (σ(a1 ), . . . , σ(ani ))> où σ : IR → IR est une fonction
sigmoïde, par exemple,
1
σ(s) = pour un λ > 0
1 + e−λs
qui a l’avantage d’être dérivable partout. Un autre exemple est la fonction ReLU (“Rectified Linear
Unit") :
σ(s) = max(s, 0)
105
On note ensuite que
∂C
∆k = .dσk , ∆j−1 = ∆j Aj dσj−1 , j = k, k − 1, . . . , 2 (5.5)
∂ak
autrement dit, à partir de ∆k , on calcule par indices décroissants (“backward phase"), ∆k−1 =
∆k Ak dσk−1 , puis ∆k−2 = ∆k−1 Ak−1 dσk−2 , etc, jusqu’à ∆1 .
Ainsi, pour calculer la différentielle de F , on commence d’abord par la phase forward, i.e., le
calcul (5.4), puis on implémente la phase backward, i.e., le calcul (5.5).
La phase backward s’appelle rétropropagation des gradients (“gradient backpropagation" en an-
∂C
glais) pour la raison suivante. La fonction C est à valeur réelles, donc la différentielle ∂a k
est
identifiée à une matrice ligne. Sa transposée est le gradient :
∂C >
= ∇ak C
∂ak
On peut écrire (5.5) sous la forme transposée :
∆> >
k = (dσk ) ∇ak C, ∆> > > >
j−1 = (dσj−1 ) Aj ∆j , j = 2, . . . , k
ce qui revient à propager en arrière le gradient. Cela est très relié aux techniques utilisées en
différentiation automatique évoquées en section 5.1.
Bien entendu, les théories de machine learning, deep learning comportent de nombreuses consi-
dérations. Elles sont actuellement en pleine évolution.
Utilisation de TensorFlow. Un peu comme AMPL est une boîte noire (extrêmement efficace) pour
coder facilement et rapidement des problèmes d’optimisation, même très difficiles, TensorFlow est
un outil libre disponible sur le web pour créer des modèles d’Intelligence Artificielle (IA) et de les
entraîner (Machine Learning, Deep Learning). Il faut coder en Python.
Pour installer TensorFlow : [Link]
On trouve sur le web quantité de bons tutoriels pour apprendre à utiliser TensorFlow.
Il existe beaucoup d’outils pour coder des procédures d’apprentissage. Le deep learning, l’IA,
sont des applications de l’optimisation qui, bien qu’étant assez anciennes, ont actuellement un
très fort potentiel applicatif, car la puissance computationnelle des ordinateurs permet désormais
d’aborder des problèmes que l’humain a du mal à appréhender. Mais ce n’est jamais que de l’opti-
misation, et l’ordinateur ne fait pas autre chose que de calculer ce qu’on lui a demandé de calculer !
106