La méthode de Newton et ses variantes pour l’optimisation

Thierry Audibert
thierry.audibert@maths2b.fr

21 décembre 2013

Disponible sur http://www.univenligne.fr sous le nom NewtonGauss.pdf

Table des matières

1 Objectifs du chapitre 2

2 Méthode de Newton en dimension 1 2

2.1 Première approche de la méthode de Newton . . . . . . . . . . . . . . . . . . . 2
2.2 Vitesse ou ordre de convergence d’un méthode itérative . . . . . . . . . . . . . . 5
2.3 De la difficulté de mise en œuvre de la méthode de Newton . . . . . . . . . . . . 9

3 Méthode de Newton en dimensions supérieures 10

3.1 Brefs rappels de calcul différentiel . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 La méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4 Méthodes d’optimisation, généralités 16

4.1 Rechercher les zéros du gradient pour optimiser . . . . . . . . . . . . . . . . . . 16
4.2 Un préalable : minimisation selon une ligne de descente . . . . . . . . . . . . . . 16

5 Problèmes de moindres carrés 24

5.1 Exemples de problèmes de moindres carrés . . . . . . . . . . . . . . . . . . . . 24
5.2 Algorithme et matrice de Gauss-Newton . . . . . . . . . . . . . . . . . . . . . . 27
5.3 Mise en œuvre : MAPLE et formule harmonique . . . . . . . . . . . . . . . . . 29
5.4 Mise en œuvre : Scilab et formule harmonique . . . . . . . . . . . . . . . . . . . 36

6 Méthode quasi-Newton (BFGS) 37

7 Annexe : historique des méthodes de moindres carrés 38

8 Correction des exercices 39

1
1 Objectifs du chapitre
La méthode de Newton consiste à approcher une solution x∗ de l’équation f (x) = 0 pour laquelle
Df (x∗ ) est inversible en la regardant comme solution de F (x∗ ) = x∗ avec

F (x) = x − D−1 f (x) ◦ f (x).

Nous introduisons la méthode en commençant par un étude détaillée (convergence, vitesse de

convergence, bassin d’attraction...) pour les fonctions d’une variable réelle. Cela permet de propo-
ser certaines démonstrations au niveau L1. L’étude de la méthode pour des fonctions de plusieurs
variables demande un minimum de calcul différentiel et se situe au niveau L2 et elle est traitée
séparément bien que les résultats soient identiques.
Plus que de résoudre des systèmes d’équations, l’objectif est ici de proposer une méthode d’op-
timisation : on recherche généralement les points réalisant le minimum d’une fonction à valeurs
réelles f en recherchant les points en lesquels g = Df (ou f 0 ) s’annule 1 . La fonction auxiliaire
s’exprime alors F (x) = x − Dg(x)−1 ◦ g(x) et fait intervenir la hessienne de f : Dg = Hf (ou
g 0 = f ”).
Ce qui fait l’intérêt de la méthode c’est la rapidité de convergence des suites xn+1 = F (xn )
lorsque convergence il y a. Mais le choix d’un x0 assurant cette convergence n’est pas aisé et
le calcul de Dg(x)−1 = Hf (x)−1 peut s’avérer coûteux. Cela nous incite à introduire comme
méthodes effectives, les méthodes Gauss-Newton pour l’étude de problèmes de moindres carrés
puis les méthodes de type quasi-Newton.
Une mise en œuvre de chacune des méthodes exposées est détaillée et implémentée sous MAPLE
et/ou Scilab pour la résolution de problèmes effectifs : détermination des coefficients de marée (et
donc prévision des hauteurs de marée), recherche du centre d’un cercle dont on connait approxi-
mativement quelques points (issu d’une problèmatique industrielle). C’est encore l’occasion de
discuter des difficultés numériques lors du passage au monde réel.

2 Méthode de Newton en dimension 1

On rappelle brièvement le principe de la méthode de Newton pour la résolution approchée d’une
équation f (x) = 0 lorsque f est de classe C 2, f 0 (x∗ ) 6= 0 en x∗ tel que f (x∗ ) = 0.

2.1 Première approche de la méthode de Newton

Théorème 1
Soient f une fonction de classe C 2 sur l’intervalle ouvert I =]a, b[, c un point de I en lequel
f (c) = 0 et |f 0 (c)| > 0 et F la fonction de classe C 1 définie par
 
f (x)
F (x) = x − 0
f (x)
Alors,
1. Ce qui, rappelons le, n’est pas une condition suffisante pour que la fonction atteigne un extremum local et encore
moins global en ce point.

2
 1. La fonction F est définie sur un voisinage ouvert de c;
2. c est un point fixe super attractif de F (ie : F (c) = c et F 0 (c) = 0).
3. Il existe δ > 0 tel que
(a) l’intervalle J = [c − δ, c + δ] est contenu dans DF et stable par f ;
(b) pour tout x0 ∈ J, la relation xn+1 = F (xn ) pour tout n ∈ N, définit une suite de J N
qui vérifie en outre
 2n −1
M2 n
|xn − c| 6 |x0 − c|2 (2.1)
2m1

avec M2 = sup |f ”(x)| et m1 = inf |f 0 (x)| .

|x−c|6δ |x−c|6δ
4. Il existe un ouvert Ω contenant c tel que pour tout x0 ∈ Ω, la suite récurrente définie par
xn+1 = F (xn ) pour tout n ∈ N, converge vers c (et il existe donc un rang à partir duquel
xn ∈ J).

Démonstration (L1)
Soient f une fonction de classe C 2 sur l’intervalle ouvert I =]a, b[ et c ∈ I tel que que |f 0 (c)| > 0.
• Comme f 0 est continue, l’ensemble des points en lesquels |f 0 (x)| > 0 est ouvert 2 3 dans I et
f ”(x)f (x)
contient c; on a bien sûr F (c) = c et F 0 (x) = est nulle en x = c.
f 02 (x)
• Comme F 0 est continue (F est de classe C 1 ) et F 0 (c) = 0, il existe δ > 0 tel que

J = [c − δ, c + δ] ⊂ DF et ∀x ∈ J, |F 0 (x)| < 1.

L’intervalle J est donc stable par F. En effet, par le théorème des accroissements finis, pour tout
x ∈ J, il existe tx ∈ J (compris entre x et c) tel que

|F (x) − F (c)| = |F (x) − c| ≤ |F 0 (tx )||x − c| < |x − c| < δ

En particulier, la relation xn+1 = F (xn ) définit une suite d’éléments de J pour tout x0 ∈ J.
• Soit x0 ∈ J et (xn )n définie par xn+1 = F (xn ).
Observons que

f (xn ) + f 0 (xn )(c − xn )

 
f (xn )
xn+1 − c = xn − 0 −c=−
f (xn ) f 0 (xn )

et par la formule de Taylor-Lagrange, il existe un point txn compris entre xn et c tel que
1
f (c) = f (xn ) + f 0 (xn )(c − xn ) + f ”(txn )(c − xn )2
2
2. Rappel de Topologie niveau L2 (en 2012) : si f est continue, l’image réciproque d’un ouvert est un ouvert
3. Si cette notion n’est pas claire pour vous, sachez qu’une preuve d’un résultat analogue, de niveau L1, est détaillée
dans la démonstration du théorème 3

3
Comme f (c) = 0, on retrouve avec M2 = sup |f ”(x)| et m1 = inf |f 0 (x)| :
|x−c|6δ |x−c|6δ

 f (xn ) + f 0 (xn )(c − xn )  |f ”(txn )|

 
M2
|xn+1 − c| = 

0
= (c − xn )2 | ≤ |xn − c|2 (2.2)
f (xn )  2 2m1

m1 et M2 sont bien définis puisqu’une fonction continue sur un segment de R est bornée (et atteint
ses bornes).
• A partir de là, on démontre par récurrence sur n la relation (2.1) :
 2n −1
M2 n
|xn − c| 6 |x0 − c|2
2m1

– Elle est immédiate pour n = 0.

– Supposons (2.1) établie pour une certaine valeur de n. En substituant dans (2.2), il vient

M2
|xn+1 − c| ≤ |xn − c|2
2m1
n !2
M2 2 −1

M2 2n
≤ |x0 − c|
2m1 2m1
 n+1
M2 2 −2+1

n+1
≤ |x0 − c|2
2m1

• On observera que la dernière inégalité n’implique pas la convergence. Par contre, comme elle
est de la forme  2n+1
M2 M2
× |x0 − c|
2m1 2m1
M2
la suite converge dès lors que |x0 − c| < 1. L’ouvert Ω =]c − δ1 , c + δ1 [ avec δ1 <
2m1
 
M2
min δ, fait donc l’affaire.
2m1


4
Théorème 2 bassin d’attraction de c
Sous les hypothèses du théorème 1, l’ensemble des points x0 ∈ I pour lesquels la suite récurrente
définie par xn+1 = F (xn ) pour tout n ∈ N, converge vers c est un ouvert de I. On l’appelle
bassin d’attraction du point fixe c.

Démonstration (L2) 4
Notons Wc l’ensemble des points a ∈ I pour lesquels la suite récurrente définie par x0 = a et
xn+1 = F (xn ) pour tout n ∈ N, converge vers c.
Considérons x0 ∈ Wc . Il existe un rang Nx0 à partir duquel xn ∈ Ω =]c − δ1 , c + δ1 [ où δ1 est
défini comme dans la démonstration du théorème 1. Observons au passage que Ω est stable par F
et donc par ses itérées.

−1
Comme la fonction itérée F (Nx0 ) est continue, O = F (Nx0 ) (Ω) est un ouvert de I. Or, pour
tout y0 ∈ O, la suite (yn )n (notation supposée évidente) est bien définie et telle que yNx0 +p ∈ Ω
pour tout p ∈ N.
Par construction de Ω,
p
M2 2 −1 

2p
|yNx0 +p − c| 6 yNx0 − c → 0
2m1
et lim yn = c. Nous avons montré que y0 ∈ O ⇒ y0 ∈ Wc .
Ainsi, Wc est voisinage de chacun de ses points et c’est bien un ouvert.


2.2 Vitesse ou ordre de convergence d’un méthode itérative

Définition 1 vitesse de convergence
On dit qu’une suite (xn )n de limite ` a une convergence d’ordre p ∈ N lorsque le quotient
|xn+1 − `|
|xn+1 − `|p
admet une limite strictement positive (ce qui suppose qu’il soit défini à partir d’un certain rang).
Lorsque p = 1 on parle de convergence géométrique ou linéaire, lorsque p = 2 on parle de
convergence quadratique.

Illustrons cette notion en montrant qu’une suite récurrente (ou système dynamique) admet une
convergence linéaire, alors que la méthode de Newton est quadratique. Le but est ici de mettre
en valeur la rapidité de convergence de la méthode de Newton qui fait tout son intérêt malgré les
problèmes que l’on rencontre pour la mettre en œuvre. Remarquons au passage que la démonstration
du théorème (3) fait apparaı̂tre tous les paradigmes 5 de l’étude des suites récurrentes et c’est pour-
quoi nous l’avons détaillée.
4. pour la notion d’image réciproque d’un ouvert
5. Consulter un dictionnaire

5
Théorème 3 ordre de convergence vers un point fixe attractif
Soient f une fonction de classe C 1 sur un intervalle ouvert I ⊂ R et a ∈ I tels que
(
f (a) =a
0
0 < |f (a)| < 1

Il existe alors α > 0 tel que

1. J =]a − α, a + α[⊂ I;
2. J est stable par f ;
3. pour tout x0 ∈ J, la suite récurrente xn+1 = f (xn ) converge vers a;
4. lorsque x0 6= a, la convergence est linéaire et

|xn+1 − a|
lim = |f 0 (a)|
|xn − a|

Démonstration (L1)
• Comme |f 0 | est continue, pour toute constante k telle que |f 0 (x0 )| < k < 1, il existe un voisi-
nage de a, J =]a − α, a + α[, contenu dans I, sur lequel 0 < |f 0 (x)| < k < 1.

Preuve de ce résultat classique que l’on pourra sauter :

|f 0 | est continue en a signifie que pour tout ε > 0, il existe αε > 0 tel que pour
tout x ∈ I,
|x − a| ≤ αε ⇒ |f 0 (x)| − |f 0 (x0 )| ≤ ε1
 

Étape 1 : On considère alors k tel que |f 0 (x0 )| < k < 1 et on pose ε1 =

1
(k − |f 0 (x0 )|) . On aura, lorsque |x − a| ≤ αε1 ,
2
1
|f 0 (x)| ≤ |f 0 (x0 )| + ε1 ≤ k + |f 0 (x0 )| ≤ k < 1


2
|f 0 (x0 )|
Étape 2 : On considère de façon analogue ε2 < et lorsque |x − a| ≤
2
αε2
0 < |f 0 (x0 )| − ε2 ≤ |f 0 (x)|
Comme I est ouvert on peut considérer α ≤ min(αε1 , αε2 ) de telle sorte que
]a − α, a + α[⊂ I et sur un tel intervalle nous avons 0 < |f 0 (x)| < k < 1.


• J est stable par f. C’est une conséquence du théorème des accroissements finis : pour x ∈ J, il
existe t compris entre x et a, tel que

|f (x) − a| = |f (x) − f (x0 )| = [f 0 (t)| |x − a| ≤ k |x − a| < α.

6
On a donc bien f (x) ∈ J. En conséquence une suite récurrente telle que x0 ∈ J et xn+1 = f (xn )
est correctement définie. De plus, xn+1 − a| = |f (xn ) − f (x0 )| ≤ k|xn − a|. De cela on déduit
que |xn − a| ≤ k|x0 − a|n ce qui assure la convergence.
• Si x0 6= a, aucun terme de la suite n’est égal à a. En effet, observons (toujours avec le théorème
des accroissements finis) que pour tout n ∈ N,
xn+1 − a = f (xn ) − f (x0 ) = f 0 (tn )(xn − a).
Si pour un certain n, xn = a, par récurrence descendante, on aurait x0 = a car sur J f 0 ne
s’annule pas.
• Ainsi, lorsque x0 6= a, nous pouvons écrire que
xn+1 − a
= f 0 (tn ) ∼ f 0 (x0 )
xn − f (x0 )
d’où la conclusion
|xn+1 − a|
lim = |f 0 (x0 )|
|xn − a|


Théorème 4 ordre de convergence de la méthode de Newton

Adjoignons l’hypothèse ”f est de classe C 3 ” à celles du théorème 1 dont nous conservons les
notations. Alors, il existe α > 0 tel que
1. J =]c − α, c + α[⊂ Wc ;
2. J est stable par F ;
3. pour tout x0 appartenant à J\{c}, la suite définie par xn+1 = F (xn ) a tous ses termes
différents de c et vérifie
|xn+1 − c| |f ”(c)|
lim =
|xn − c|2 2|f 0 (c)|
4. lorsque f ”(c) 6= 0, une telle suite admet une convergence d’ordre 2 vers c.

Démonstration (L1)
• Rappelons que F est définie sur un ouvert contenant c avec
f (x)
F (x) = x − ;
f 0 x)
f (x) f ”(x)
F 0 (x) = ;
(f 0 (x))2
f ”(x)(f 0 (x))2 − 2f (x)(f ”(x))2 + f (x)f (3) (x)f 0 (x)
F ”(x) = ;
(f 0 (x))3
f ”(c)
Au point c nous retrouvons F 0 (c) = 0, F ”(c) = puisque f (c) = 0.
f 0 (c)
• Nous pouvons considérer un réel α > 0 tel que

7
– J =]c − α, c + α[⊂ Wc (car Wc est un ouvert) ;
– |F 0 (x)| < 1 pour tout x ∈ J (car F 0 est continue et F 0 (c) = 0).
• J est alors stable par F 6 et pour x0 ∈ J, la relation xn+1 = F (xn ) définit une suite d’éléments
de J qui converge vers c (J est contenu dans le bassin d’attraction de c).
• Si x0 6= c, aucun terme de la suite n’est égal à c.
En effet, observons (toujours avec le théorème des accroissements finis que
nous pouvons invoquer car J est un intervalle contrairement à Wc peut-être)
que pour tout n ∈ N,

xn+1 − c = f (xn ) − f (c) = f 0 (tn )(xn − c).

Si pour un certain n, xn = c, par récurrence descendante, on aurait x0 = c car

f 0 ne s’annule pas là où F est définie.
• Par la formule de Taylor-Lagrange, il existe tn ∈ J compris entre xn et c tel que
1
xn+1 = F (xn ) = F (c) + F 0 (c)(xn − c) + F ”(tn )(xn − c)2
2
Il vient donc
1
xn+1 − c = F ”(tn )(xn − c)2
2
xn+1 − c 1
Nous pouvons alors écrire = F ”(tn ) d’où
(xn − c)2 2
xn+1 − c 1 f ”(c)
lim = lim F ”(tn ) =
(xn − c) 2 2 2 f 0 (c)
Lorsque cette limite est non nulle, la convergence est exactement d’ordre 2.

Une illustration numérique de la vitesse de convergence
2 x3
On considère ici la fonction f : x → x3 − x associée à la fonction de Newton F : x → .
3 x2 − 1
xn − c
Chaque tableau présente les valeurs x0 , ..., x6 ainsi que les quotients . Lorsque le
(xn−1 − c)2
f ”(c) −3 f ”(c)
point fixe est c = 1, 0 = et pour c = 0, 0 = 0.
2f (c) 2 2f (c)

xn − c xn − c
n xn n xn
(xn−1 − c)2 (xn−1 − c)2

0 2.0 * 0 0.3333333333 *
1 1.454545455 0.4545454545 1 - 0.1111111111 - 1.0
2 1.151046789 0.7310664606 2 0.002849002849 0.2307692308
3 1.025325929 1.110049608 3 - 0. 04625079707 - 0.005698144450
4 1.000908452 1.416351963 4 1.978735114 0−22 0. 09250159415
5 1.000001235 1.496827212 5 −1.549504984 × 10−65 −3.957470228 × 10−22
6 1.0 1.499995676 6 7.440616629 × 10−195 3.099009969 × 10−65

6. Toujours le théorème des accroissements finis comme dans les démonstrations des théorèmes 1 et 3

8
2.3 De la difficulté de mise en œuvre de la méthode de Newton
Comme on le voit avec l’exemple qui suit, rien ne garantit qu’une suite (xn )n donnée par la
méthode de Newton et de premier terme x0 converge.

x3
F IGURE 1 – f : x → x3 − x F : x → 2
3 x2 − 1

On trouvera de jolies études détaillées sur les bassins d’attraction pour la méthode de Newton mise
en œuvre avec des fonctions élémentaires dans l’ouvrage de Holmgren [6].
A FINIR – A FINIR – A FINIR – A FINIR –

9
3 Méthode de Newton en dimensions supérieures
Nous nous intéressons ici à la résolution de systèmes de n équations à n inconnues ou, ce qui
revient au même, à la recherche de zéros de fonctions de Rn dans lui-même. Les résultats fonda-
mentaux sont les mêmes que pour le cas des fonctions d’une variable réelle.

3.1 Brefs rappels de calcul différentiel

Nous rappelons quelques notions de calcul différentiel et quelques notations qui nous seront utiles.

Dans ce qui suit Ω est un ouvert de Rn et f une application de Ω dans Rm .

Fonctions de classe C k : soit k ∈ N∗ ; on dit que f est de classe C k si toutes ses dérivées par-
tielles d’ordre 1 à k sont définies et continues sur Ω.
Matrice jacobienne : lorsque f est de classe C 1 sa matrice jacobienne en a ∈ Ω est la matrice
DF (a) ∈ Mm,n définie par
∂f1 (a) ∂f1 (a) ∂f1 (a)
 
 ∂x1 . . .
∂x2 ∂xn 
 ∂f2 (a) ∂f2 (a) ∂f2 (a) 
 
 ... 
Df (a) =  ∂x. 1 ∂x2
..
∂xn 
.. 
 .. . . 
 
 ∂fm (a) ∂fm (a) ∂fm (a) 
...
∂x1 ∂x2 ∂xn
L’application linéaire de Rn dans Rm qui lui est canoniquement associée est la différentielle
de F en a.
Théorème de Schwarz : lorsque f est de classe C 2 , pour tous couples d’indices (i, j) ∈ [1, n]2 ,
∂ 2 f (x1 , ..., xn ) ∂ 2 f (x1 , ..., xn )
=
∂xi ∂xj ∂xj ∂xi

Hessienne : lorsque f est de classe C 2 à valeurs réelles, sa matrice hessienne en a ∈ Ω est la

matrice définie par
 2
∂ f (a) ∂ 2 f (a) ∂ 2 f (a)

 ∂x1 ∂x1 ∂x1 ∂x2 ...
 2 ∂x1 ∂xn 
 ∂ f (a) ∂ 2 f (a) ∂ 2 f (a) 

 ... 
 ∂x2.∂x1 ∂x2.∂x2
Hf (a) =  ∂x2 ∂xn 
.. 

 .. .. . 

 ∂ 2 f (a) ∂ 2 f (a) ∂ 2 f (a) 
...
∂xn ∂x1 ∂xn ∂x2 ∂xn ∂xn
Nous verrons, avec la formule de Taylor à l’ordre 2 qu’il est plus parlant de considérer cette
matrice symétrique (par le théorème de Schwarz) comme associée à une forme quadratique.
Différentielle seconde : supposons encore f de classe C 2 à valeurs réelles, son gradient g : Ω →
Rn (pour le produit scalaire canonique) est de classe C 1. On vérifie sans peine que la matrice
jacobienne du gradient dans la base canonique est aussi la hessienne de f (Dg = Hf ).
On prendra garde aux changements de repères, si toutefois on est amené à en envisager.

10
Développements limités et formules de Taylor :
• une fonction de classe C 1 f : Ω ⊂ Rn → f (x) ∈ Rm admet un DL1 en tout point a ∈ Ω
donné par :

Pp ∂
f (a + h) =h→0 f (a) + i=1 f (a)hi + o (| |h||) =h→0 f (a) + Df (a).h + o (| |h||)
∂xi

• Une fonction de classe C 2 à valeurs réelles, f : Ω ⊂ Rn → f (x) ∈ R admet un DL2 en

tout point a ∈ Ω donné par :

1t
f (a + h) =h→0 f (a) + Df (a).h + hHf (a)h + o(||h||2 )
2

• Pour une fonction de classe C 2 à valeurs dans Rm , f : Ω ⊂ Rn → f (x) ∈ Rm , on

réécrira la formule composante par composante ce qui donne :
t   
hHf1 (a)h o1 (1)
1 ..  . 
f (a + h) =h→0 f1 (a) + Df (a).h +   + ||h||2  .. 

2 t .
hHfm (a)h om (1)

1
On notera aussi f (a + h) =h→0 f (a) + Df (a).h + Df2 (a)(h, h) + o(||h||2 )
2

Changements de repères :
Rappelons les formules de changements de bases puis de repère pour les points d’un espace
affine, pour les applications linéaires et les formes quadratiques.
Soient E un espace affine, R = (O, (i, j, k)) = (O, B) et R0 = (A, (u, v, w)) = (A, B 0 )
deux repères de E. Pour M ∈ E on note XM et XM 0 éléments de Rn , coordonnées
(−les
−→
OM = x~i + y~j + z~k
respectives de M dans ces repères, ce qui signifie que −−→
AM = x0 ~u + y 0~v + z 0 w~
• On note P ∈ GLn (R) la matrice de passage de la base B à B 0 . On sait que si X et X 0
sont les coordonnées d’un même vecteur de E ~ dans ces bases on a P X 0 = X (formule de
changement de base pour un vecteur).
−−→ −→ −−→
• De la relation vectorielle AM = AO + OM , on exprime les coordonnées de ce vecteur
sans la base B, il vient (formule de changement de repère pour un point) :

0 = −X + X
P XM 0 +X =X
ou encore P XM
A M A M

• Supposons l’espace affine normé et considérons une fonction scalaire F : U ⊂ E → R,

définie sur un ouvert U ⊂ E. On associe à cette fonction F la fonction f : Ω ⊂ Rn → R
(n = dim E) définie par f (XM ) = F (M ) pour tout M ∈ E; f est donc l’expression de
F dans le repère (0, B).
En notant g(XM 0 ) = F (M ) l’expression de F dans l’autre repère, la relation qui suit entre

11
 f et g donne une formule de changement de repère pour les applications scalaires :
g(XM 0 ) = f (P X 0 + X ) = f (X ) ou encore
M A M

g(X 0 ) = f ◦ Φ(X 0 ) u, avec Φ(X 0 ) = P X 0 + XA

• Ecrivons alors les formules de Taylor pour f et g à l’ordre 2 (en supposant F, donc f et g
de classe C 2 ) :
1t
f (x0 + h) = f (x0 )) + Df (x0 ).h + hHf (x0 )h + o(||h||2 ) (3.1)
h→0 2
1t
g(x00 + h0 ) 0= g(x00 ) + Dg(x00 )).h0 + h0 Hg (x00 )h0 + o(||h0 ||2 ) (3.2)
h →0 2
Bien évidemment, avec g(x00 + h0 ) = f ◦ φ(x00 + h0 ) = f (P (x00 + h0 ) + XA ) en remplaçant
x0 = φ(x00 ) par P x00 + XA et h par P h0 (changement de base pour un vecteur), on obtient

1 t 0 t
f (x0 + h) 0= f (φ(x00 )) + Df (φ(x00 )).P h0 + h ( P Hf (φ(x00 ))P ) h0 + o(||P h0 ||2 )
h →0 2
(3.3)
1t 0
= g(x00 ) + Dg(x00 ).h0 + h Hg (x00 )h0 + o(||h0 ||2 ) (3.4)
0
h →0 2
En identifiant les termes de même ordre (unicité du DL par exemple) on obtient :

Dg(y) = Df (φ(y)) ◦ P

Hg (y) = tP ◦ Hf (φ(y)) ◦ P

Nous aurons besoin de certains des résultats suivants pour la démonstration du théorème 6. Il
supposent une bonne compréhension de la notion de différentielle (et ce n’est pas vraiment du
L2).

Théorème 5
Soit J : M ∈ GLn (R) → M −1 ∈ GLn (R).
1. J est une application de classe C ∞ de l’ouvert GLn (R) dans lui-même ;
2. ** pour toute matrice A ∈ GLn (R), la différentielle de J au point A est l’application
linéaire dJ (A) ∈ L (Mn (R)) définie par :

dJ (A)(H) = A−1 × H × A−1 (3.5)

Démonstration proposée dans l’exercice 1

12
Exercice 1 démonstration du théorème 5
1. Soient I un intervalle de R, A : I → Mn (R) et B : I → Mn (R) deux application de
classe C 1 . Montrer que l’on a

d
(A(t) × B(t)) = A0 (t) × B(t) + A(t) × B(t)
dt
Qu’en est il si B : I → Rn ?
2. On note J l’application qui à une matrice inversible associe son inverse J : M ∈ GLn (R) →
M −1 ∈ GLn (R).
(a) Montrer que GLn (R) est un ouvert de Mn (R).
(b) Justifier que J est une application de classe C ∞ de l’ouvert GLn (R) dans lui-même.
(c) Les idées claires : de quoi parlons nous ?
La différentielle de J au point A est une application linéaire : quels sont les espaces
de départ et d’arrivées ce cette application ? Quels sont les coefficients de sa matrice
jacobienne exprimées une base canonique bien choisie ? Combien y en a-t-il ?
Rappeler la définition des dérivées partielles en un point A ∈ GLn (R) dans cette base.
De quel type d’objet s’agit il ?
3. Pour nous faciliter la vie, nous allons plutôt déduire l’expression de dJ (A)H de l’étude
d’un courbe paramétrée à valeurs dans GLn (R).
(a) Soit Γ une courbe paramétrée de classe C 1 de I, intervalle de R, à valeurs dans
d
GLn (R). Exprimer J ◦ Γ(t).
dt
(b) En considérant une courbe paramétrée t → Γ(t) × Γ(t)−1 exprimer dJ (Γ(t))Γ0 (t) et
en déduire dJ (A)H lorsque A est une matrice inversible quelconque, H une matrice
quelconque.
voir correction en (8.0.1)

Exercice 2 un petit calcul préliminaire

Nous avons établi dans l’exercice 1 une formule de dérivation de t → A(t) × B(t) lorsque
A : I → Mn (R) et B : I → Mn (R) sont deux courbes paramétrées de classe C 1 .
On se propose maintenant de déterminer la différentielle d’une application x ∈ Ω → A(x)×V (x)
lorsque lorsque A : Ω → Mn (R) et V : Ω → Rn sont de classe C 1 sur Ω ouvert de Rn .
1. Dériver l’application t → A(X(t)).V (X(t)) lorsque X est une courbe paramétrée de
classe C 1 à valeurs dans Ω. Préciser avec soin un ensemble d’appartenance de chaque ob-
jet/expression introduite par le calcul.
2. En déduire la formule d(A.V )(x) = dA(x).V (x) + A(x)DV (x)
voir correction en (8.0.2)

13
3.2 La méthode de Newton
Dans le théorème qui suit, en tout point analogue au théorème 1 par sa forme, on désigne par || ||
une norme de Rn et par ||| ||| la norme matricielle subordonnée qui lui est attachée. On rappelle
que pour A ∈ Mn (R),
||AX||
|||A||| = sup
X6=0 ||X||

Théorème 6
Soient Ω un ouvert de Rn , f : Ω → Rn , une fonction de classe C 2 , c un point de Ω en lequel
f (c) = 0 et Df (c) est inversible. On considère F la fonction de classe C 1 définie par

F (x) = x − [Df (x)]−1 .f (x)

Alors,
1. La fonction F est définie sur un voisinage ouvert de c;
2. c est un point fixe super attractif de F (ie : F (c) = c et Df (c) = 0).
3. Il existe δ > 0 tel que
(a) la boule fermée B = B̄(c, δ) est stable par F ;
(b) pour tout x0 ∈ B, la relation xn+1 = F (xn ) pour tout n ∈ N, définit une suite de B N
qui vérifie en outre
n −1 n
||xn − c|| 6 K 2 ||x0 − c||2 (3.6)

pour un certain K > 0.

4. Il existe un ouvert Ω contenant c tel que pour tout x0 ∈ Ω, la suite récurrente définie par
xn+1 = F (xn ) pour tout n ∈ N, converge vers c (et il existe donc un rang à partir duquel
xn ∈ B).

Remarques
1. On peut préciser ce résultat en montrant que l’on peut choisir pour tout ε > 0, B et K de
telle sorte que
 
1 −1 2
K < max |||Df (x)] ||| × max ||D f (x)(h, h)|| + ε
2 x∈B x∈B,||h||=1

2. Les théorèmes 2 et 4 concernant le bassin d’attraction et la vitesse de convergence s’énoncent

et se démontrent de façon analogue en dimension quelconque
Démonstration (L2)
La démonstration que nous proposons ici en exercice, qui peut être sautée, fait largement appel
au calcul différentiel avec les résultats des exercices 1et 2. Elle a néanmoins l’avantage d’être
proche de la preuve en dimension 1. On peut en trouver une autre, qui repose sur un calcul de
développement limité dans l’ouvrage de Demailly [5] par exemple.

14
Exercice 3 démonstration du théorème 6
Soit f comme dans le théorème et F (x) = x − [Df (x)]−1 .f (x).
1. Justifier que F est définie sur un ouvert contenu dans Ω, calculer DF (x) en réécrivant

F (x) = x − (J ◦ df )(x).f (x)

Vérifier que F (c) = c et que dF (c) = 0.

2. En déduire l’existence d’une boule fermée B = B̄(c, δ) , contenue dans le domaine de
définition de F telle que
R1 d
(a) B est stable par F (penser à écrire F (c + h) − F (c) = 0 F (c + t h) dt);
dt
(b) Pour tout x0 ∈ B, la suite de B N définie par xn+1 = F (xn ) pour tout n ∈ N, vérifie
||xn+1 − c|| ≤ K||xn − c||2 avec K que l’on comparera à une expression proche de

max |||Df (x)]−1 ||| × max ||D2 f (x)(h, h)||

x∈B x∈B,||h||=1

indication : penser à utiliser la formule de Taylor

voir correction en 8.0.3

15
4 Méthodes d’optimisation, généralités
4.1 Rechercher les zéros du gradient pour optimiser
Soit φ : Ω ⊂ Rn → R, une fonction numérique de classe C 3 définie sur un ouvert de Rn . On sait
−−→
que si φ atteint en a ∈ Ω un minimum local ou global, grad(f )(a) = 0.
−−→
La fonction g = grad(f ) est une application de Ω dans Rn de classe C 2 et il peut paraitre
avantageux de rechercher les points en lesquels g s’annule avec la méthode de Newton mais cela
pose plusieurs problèmes :
−−→
– grad(f )(a) = 0 est une condition nécessaire, non suffisante pour que f admette un minimum
local en a et il faudra pour chaque point stationnaire ainsi déterminé s’assurer qu’on y rencontre
bien un minimum local (à moins qu’une étude préalable n’ait déjà clarifié le problème) ;
– à cela s’ajoutent les difficultés de mise en œuvre de la méthode elle-même :
- la convergence n’est pas garantie et dépend du choix de x0 ;
- la méthode revient à itérer la fonction

G(x) = x − [Dg(x)]−1 g(x) = x − [Hf (x)]−1 g(x) (4.1)

où Dg, jacobienne du gradient de f, est en fait la matrice hessienne de f ; il faut donc calculer
g, calculer Hf et ensuite l’inverser. On comprend que la difficulté n’est pas la même selon que
l’on sait calculer g et Dg = Hf formellement ou pas.

On est alors conduit à adapter cette méthode en remplaçant l’itération xn+1 = G(xn ), qui se récrit
xn+1 = xn − [Hf (xn )]−1 g(xn ), à partir des idées suivantes :

• on remplace [Hf (xn )]−1 par une matrice Wn définie positive dont le calcul et l’inversion sont
plus économiques ;
• on cherche un réel λ∗ qui réalise un minimum local de t → f (xn − tWn g(xn )) sur [0, 1] et on
pose

xn+1 = xn λ∗ Wn g(xn ) (4.2)

ce qui nous garantit que la suite (f (xn )n décroit et nous permet d’espérer une convergence qui
n’était pas gagnée.

On parle alors de méthodes quasi-newton. Nous donnons une méthode de recherche selon une
ligne de descente en (4.2) et présentons deux méthodes quasi-newton en détails en (5) et en (6)
avec une mise en œuvre de chacune d’elles sur un problème effectif.

4.2 Un préalable : minimisation selon une ligne de descente

Nous serons amenés, pour minimiser une fonction de plusieurs variables f : Ω → R, à définir des
suites approximantes (xn )n pour lesquelles le calcul de xn+1 , xn étant donné, sera défini en trois
étapes :
1. On recherche une direction de descente dn (ie : telle que λ → f (xn + λdn ) & sur un
voisinage à droite de 0).

16
 −−→
• Une première idée serait de choisir par exemple dn = −grad(f )(xn ), ce qui se comprend
−−→
bien en écrivant un DL1 , f (xn + h) = f (xn )+ < grad(f )(xn )|h > +o(h), où l’on voit
que le second terme est minimal pour un module ||h|| constant lorsque l’angle entre h et le
gradient est égal à π...
−−→
• Nous verrons aussi l’intérêt de poser dn = G−1 grad(f )(xn ) avec H symétrique positive...
2. Comme rien ne garantit que le point xn+1 = xn +dn réalise f (xn+1 ) ≤ f (xn ), on recherche
un scalaire λ∗ appartenant à un certain intervalle (= [0, +∞[) tel que

f (xn + λ∗ dn ) = min f (xn + λdn )

3. On pose alors xn+1 = xn − λ∗ dn .

Avant de préciser de telles méthodes, intéressons nous au deuxième point et montrons comment
on peut rechercher un minimum local en dimension 1. Les exercices 4 et 5 constituent un bon
échauffement. On propose ensuite une description de la méthode de Wolfe.

Exercice 4 un préliminaire à but didactique (L1)

Soit f : I → R une fonction définie sur un intervalle I de R.
1. Montrer que f est convexe sur I si et seulement si pour tous x1 < x2 ≤ x3 < x4 dans I,
f (x2 ) − f (x1 ) f (x4 ) − f (x3 )
≤
x2 − x1 x4 − x3

2. Montrer que si f est strictement convexe et continue sur I = [a, b], elle y admet un minimum
atteint en un seul point. Que devient l’inégalité précédente lorsque les points x1 , ..., x4 sont
distincts ?
3. On suppose maintenant que f est strictement convexe et continue sur I = [a, b] et on note
µ le point en lequel f atteint son minimum.
On se donne t ∈]0, 1/2[ et on partage I = [a, b] en trois intervalles de la façon suivante :
a+b


x1 = − t (b − a)

 2
I = [a, x1 ] ∪ [x1 , x2 ] ∪ [x2 , b] avec
x2 = a + b + t (b − a)



2
(a) On note p0 , p2 , p2 les pentes de f sur ces trois intervalles. Montrer les implications
suivantes :
– p0 < p1 < p2 < 0 =⇒ µ ∈ [x2 , b];
– p0 < p1 < 0 < p2 =⇒ µ ∈ [x1 , b];
– p0 < 0 < p1 < p2 =⇒ µ ∈ [a, x2 ];
– 0 < p0 < p2 < p2 =⇒ µ ∈ [a, x1 ];
(b) Déduire de cela un algorithme de calcul approché du minimum de f (programmez le
effectivement Maple, Python, Scilab...).
(c) Réduire et estimer la complexité de l’algorithme proposé (nombre d’itérations ou
nombre de calculs de f en un point pour une approximation à 10−n près par exemple).
voir corrigé en 8.0.4

17
Exercice 5 Recherche de minimums, règle d’or et gradient (L2) 7
1. Soit I = [a, b] un intervalle de R et t0 ∈ I. On considère dans cette question une fonction
f : I → R telle que

f & strictement sur ] − ∞, t0 ] ∩ I et f % strictement sur [t0 , +∞[∩I.

(a) On définit deux suites (an )n et (bn )n de la façon suivante :

– [a0 , b0 ] = [a, b];
– On suppose an et bn déjà définis et on pose λn = an + r2 `n , µn = an + r`n avec r
solution positive de x2 + x = 1 et `n = bn − an .
Si f (λn ) ≥ f (µn ) alors an+1 = λn et bn+1 = bn
Si f (λn ) ≤ f (µn ) alors an+1 = an et bn+1 = µn
Démontrer que (an )n et (bn )n convergent vers le point en lequel f atteint son mini-
mum. Majorer `n .
(b) Programmer cet algorithme
 et le tester avec f (x) = sin x sur I = [π, 2 π] puis avec la
3 − x si 0 ≤ x ≤ 1

fonction g(x) = x + 1 si 1 ≤ x ≤ 2 sur I = [0, 3]. Dans chaque cas, donner un

3 si 2 < x ≤ 3

encadrement de longueur 10−7 du réel en lequel le minimum sur I est atteint.
2. On s’intéresse maintenant à la recherche de minimum de fonctions de plusieurs variables :
on suppose que F est définie sur R2 , à valeurs réelles, de classe C 1 et que sa restriction à
une droite quelconque est minorée.
(a) Soient X ∈ R2 , et U ∈ R2 , un vecteur unitaire. On pose φ(λ) = F (X+λU ); exprimer
φ0 (0) en fonction des vecteurs U et grad(F )(X). Pour quel choix de U, φ0 (0) est elle
minimale ?
(b) Pour X0 ∈ R2 , on définit une suite de points (Xn )n par récurrence en posant :
• si grad(F )(Xn ) 6= 0,
−grad(F )(Xn )
– Un =
||grad(f )(Xn )||
– Xn+1 = Xn + λ∗n Un où λ∗n réalise le minimum de F (Xn + λUn );
• si grad(F )(Xn ) = 0, Xn est le dernier terme de la suite qui est donc finie.
Montrer que si la suite (F (Xn ))n est définie sur N, elle est convergente.

(c) On donne à titre d’exemple F (x, y) = ln 1 + x2 + y 2 + xy . Calculer son gradient ;

écrire un programme permettant le calcul des 12 premiers termes de la suite (xn , yn )

(on pourra faire appel à la fonction de la première question).
voir corrigé en 8.0.5

7. Posé à peu près sous cette forme à l’Oral Centrale-Supélec, épreuve avec Maple ou Mathematica— 2011

18
Procédé de Wolfe

Soit q : t → q(t) une fonction définie sur [0, +∞[ telle que q 0 (0) < 0,dont on souhaite évaluer
grossièrement 8 un minimum local (en pratique ce sera q(t) = f (xn + t dn )).
On définit pour cela le processus :
– on se donne m1 , m2 tels que 0 < m1 < m2 < 1;
– on se donne une valeur t > 0;
– on itère la règle de Wolfe ci-dessous, en choisissant à chaque étape un t > 0 compris entre tmin
et tmax (avec la convention que si tmax = 0, tmin < t < +∞)

q(t) − q(0)
si > m1 q 0 (0) alors tmax := t (on cherche plus petit)


t − 0








q(t) − q(0)

si ≤ m1 q 0 (0) et q 0 (t) ≥ m2 q 0 (0) on s’arrête là


 t − 0



q(t) − q(0)


≤ m1 q 0 (0) et q 0 (t) < m2 q 0 (0) alors tmin := t (on cherche plus grand)

 si

t−0
8. Grossièrement, car c’est le minimum de f qui nous intéresse : c’est le procédé global qui doit nous permettre de
nous en rapprocher alors qu’un calcul fin du minimum de q à chaque itération serait inutilement pénalisant.

19
Illustration graphique de la règle de Wolfe
Les arguments, q et q1 désignent respectivement la fonction à minimiser et sa dérivée, t est
la valeur à tester...

Wolfe:=proc(q,q1,t,m1,m2)
local tv, d0,dt;

tv :=(q(t)-q(0))/((t-0));
d0 :=q1(0);
dt :=q1(t);

if evalf(tv-m1*d0)>0
then ga #pour aller à gauche
elif evalf(dt-m2*d0)>=0 then fin
else droite #pour aller à droite
fi;
end:

q := x-> 1-6*x*sin(3*x+1/2)/(1+xˆ2);
q1 := unapply(diff(q(t), t), t);
plot({m1*q1(0)*t+q(0), m2*q1(0)*t+q(0), q(t)}, t=0..4,-4..4);

Matrix([[seq(4*t/(10.), t = 1 .. 10)],
[seq(Wolfe(q, q1, 4*t*(1/10), .3, .8), t = 1 .. 10)]])

 
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4...
dr dr fin fin ga ga ga ga ga ga ga fin ga ga ga ga ga...

20
• Algorithme : avec la procédure Wolfe du tableau précédent, une forme possible de l’algo-
rithme (les variantes dépendent du domaine de définition de q (est il borné ou pas ?), de la façon
d’initialiser t, d’extrapoler lorsque tmax = 0, d’interpoler, etc...

LineSearchWolfe := proc (q, q1, t0, m1, m2)

local arret, tmax, tmin, t, w;

t := t0; # pour la preuve t0>0

arret := false;
tmin := 0;
tmax := 0;

while arret = false do

w := Wolfe(q, q1, t, m1, m2);
if w = fin
then
arret := true
elif w = dr
then tmin := t;
if tmax = 0
then t := 10*t
else t := (1/2)*t+(1/2)*tmax
# ou evalf((1/2)*tmin+(1/2)*tmax)
end if
else tmax := t;
t := (1/2)*tmin+(1/2)*tmax
# ou evalf((1/2)*tmin+(1/2)*tmax)
end if
end do;
t
end proc

21
 Une trace avec la même fonction q que pour le graphe... Il va de soi que pour des exemples effectifs
on travaille en flottants et non en calcul formel comme ici.

LineSearchWolfe(q, q1, 1, m1, m2);

{--> enter LineSearchWolfe, args = q, q1, .4, .9
5
false
0
0
ga
5
5/2
ga
5/2
5/4
ga
5/4
5/8
fin
true
5/8
<-- exit LineSearchWolfe (now at top level) = 5/8}
5/8

Théorème 7 terminaison
Soit q une fonction définie et dérivable sur [0, +∞[ telle que
– q 0 (0) < 0;
– inf q(t) > −∞.
Alors, l’algorithme décrit par la fonction LineSearchWolfe termine.

Démonstration (L1+)
Si l’algorithme ne termine pas, il y a une infinité d’appels à la fonction Wolfe qui retournent soit
’dr’ soit ’ga’.
• Supposons qu’il y ait une infinité d’appels retournant ’dr’ sans que tmax ne soit jamais affecté.
q(t) − q(0)
Alors, t(k) → +∞ et comme ≤ m1 q 0 (0) < 0, cela impose lim q(t(k) = −∞ ce qui
t−0
est en contradiction avec les hypothèses.

22
• Dans les autres cas, à partir d’un certain rang, l’exposant (k) signifiant que la variable est évaluéé
(k) (k)
à la fin de la k ième itération, on a 0 ≤ tmin < t(k) < tmax où les suites encadrantes sont adjacentes
(k) (k)
(par construction). Notons t∗ = lim tmin = lim t(k) = lim tmax .
Comme les affectations à tmax ou tmin proviennent des deux règles

q(t) − q(0)
si > m1 q 0 (0) alors tmax := t (Wolfe retourne ’ga’)


t−0




q(t) − q(0)


≤ m1 q 0 (0) et q 0 (t) < m2 q 0 (0)

 si
 alors tmin := t (Wolfe retourne ’dr’)
t−0
les relations qui suivent sont invariantes à partir de ce rang
0
q(t(k) (k)
max ) > q(0) + m1 tmax q (0) (4.3)
(k) (k) (k)
q(tmin ) ≤ q(0) + m1 tmin q 0 (0) et q 0
(tmin ) 0
< m2 q (0) (4.4)

ce qui, par passage à la limite, donne

q(t∗ ) = q(0) + m1 t∗ q 0 (0) et q 0 (t∗ ) ≤ m2 q 0 (0) (4.5)

(k) (k)
q(tmax ) − q(t∗ ) q(tmin ) − q(t∗ )
Cela nous permet de comparer les taux de variations (k)
et (k)
à m1 q 0 (0).
tmax − t∗ −
tmin t∗
Il vient alors à la fois q 0 (t∗ ) = m1 q 0 (0) et q 0 (t∗ ) ≤ m2 q 0 (0) ce qui est contradictoire puisque
q 0 (0) < 0.

Que pouvons nous dire de q au point t ainsi obtenu ? Voyons cela lorsque q est de la forme q(t) =
f (xn + t dn ) avec des hypothèse de régularité sur f et sur son gradient, ainsi que sur le choix de
(dn )n .

23
5 Problèmes de moindres carrés
5.1 Exemples de problèmes de moindres carrés
Un problème de moindres carrés est un problème de la forme :
p p
1X 2 1X
minimiser f (x) avec f (x) = rk (x) = (φ(x, tk ) − zk )2
2 2
k=1 k=1

où chaque fonction rk2 (x) = (φ(x, tk )−zk )2 représente le carré de l’écart entre la valeur théorique
d’une fonction φ dépendant des paramètres x au point tk et la valeur zk mesurée en ce même tk .
On cherche généralement dans un tel cas à estimer les paramètres inconnus à partir d’un nombre
plus ou moins grand de relevés expérimentaux. Illustrons cela de quelques exemples.

Déterminer les coefficients de marée

La hauteur de la marée dans un lieu donné (un port par exemple) dépend des caractéristiques
astronomiques que sont les trajectoires respectives de la Terre, de la Lune et du Soleil ainsi
que de la configuration des fonds marins et de la forme des côtes 9 .
Cette hauteur est donnée par une fonction de la forme :
N
X
h(t) = φ(a1 , ..., aN , ν1 , ..., νN , t) = h0 + ai cos (ωi t + νi )
i=1

Cette expression (formule harmonique), dès lors que l’on dispose de relevés en nombre
suffisant, ce qui est le cas, se prête bien à un étude numérique. Elle est exploitée par le
service hydrographique et océanographique de la marine depuis 1992.
Les pulsations ωi sont universelles car ce sont celles des potentiels des forces de marées
d’origine lunaire et solaire. Elles sont reconnues par analyse de Fourier de relevés sur de
nombreuses années. h0 se mesure facilement car c’est la hauteur moyenne de l’eau, les
paramètres (ai )i et (νi )i dépendent du lieu considéré et sont les inconnues que l’on veut
évaluer.
9. bien entendu, la hauteur de l’eau dépend aussi de phénomènes météos : température, vents, pression at-
mosphérique...

24
 Pour déterminer les paramètres dont la connaissance permet de prévoir la hauteur d’eau à
tout instant, on recherche ceux qui réalisent le minimum de la fonction de 2N variables :
p p
1X 1X 2
f (a1 , ..., aN , ν1 , ..., νN ) = (φ(a1 , ..., νN , tk ) − zk )2 = rk (a1 , ..., νN )2
2 2
k=1 k=1

où les (zk )1≤k≤p sont les hauteurs observées en des instants (tk )1≤k≤p .
Pour les sites où la marée est importante, la prise en compte de plus d’une centaine d’ondes
composantes (ai cos(ωi t + νi )) peut être nécessaire pour un calcul précis. Ainsi la fonction
à minimiser f, peut être une fonction de 2N ≈ 200 variables.
Les paramètres sont déterminés en analysant les observations réalisées pendant une longue
période (au minimum 29 jours, correspondant à une lunaison et, là où la forme de la marée
est particulièrement complexe, une année). On a alors au moins de l’ordre de P ≈ 700
termes dans la somme.
Références :
– le site du SHOM [7] donne quelques brèves explications sur la formule harmonique (re-
prises dans Wikipedia).
– l’ouvrage très complet [8] de Bernard Simon (Éditions de l’Institut océanographique
de Paris) aborde les aspects physiques de l’étude de la marée, explique comment on
détermine les pulsations ωi ...

Estimer le centre d’un cercle

On veut reproduire un objet industriel de bord circulaire pour lequel on a pu relever les côtes
de quelques points situés sur la circonférence (figure). Bien entendu, l’objet est abimé, les
mesures entachées d’erreurs mêmes minimes. Comment déterminer le rayon au mieux, où
placer le centre connaissant approximativement ces points ?

Là encore on se ramène à un problème de moindres carrés :

– on note φ(A, r, M ) la projection du point M sur le cercle C(A, r);
– en notant M1 , M2 , ...Mp les p points dont les cotes ont été relevées, on considère la
fonction d’erreur
p
1X
f (A, r) = ||φ(A, r, Mk ) − Mk ||2
2
k=1

25
 On recherche alors les paramètres (xA , yA , r) qui réalisent le minimum de cette fonction.

Paramétrage d’un potentiel empirique

Un exemple en économie

Les moindres carrés linéaires

A la suite d’un relevé expérimental, on dispose de N données (xk , yk ) ∈ R2 et on recherche
une fonction f, combinaison linéaire de d + 1 fonctions données (ui )0≤i≤d (ie : f (x) =
Pd t
i=0 ai ui (x)) telle que le vecteur F = [f (x1 ), f (x2 ), ..., f (xN )] réalise le minimum de

N
X
||Y − F ||22 = |yk − f (xk )|2
k=1

avec Y =t |y1 , y2 , ..., yN ].

La solution de ce problème est proposée dans l’exercice qui suit.
Exercice 6 moindres carrés linéaires
1. (a) Expliciter un sous-ev V de RN auquel appartient le vecteur

F =t [f (x1 ), f (x2 ), ..., f (xN )]

lorsque f est combinaison linéaire des (ui )i . Montrer que l’on peut l’écrire F
sous la forme F = U A où U est une matrice que l’on décrira avec soin, A =t
[a0 , a1 , ..., ad ]. Exprimer la dimension de V en fonction de la matrice U.
(b) Approximation polynomiale : on choisit pour les fonctions ui , ui (x) = xi . Quelle
est la matrice U dans ce cas ?
2. Justifier que l’existence d’un vecteur F0 qui réalise le minimum de ||Y −F || est assurée
quelque soit la norme || || choisie dans RN . Pour la suite du problème, on choisit la
norme euclidienne dans RN définie par
v
uN
uX
||X||2 = (t X.X)1/2 = t x2i
i=1

3. On pose J(A) = ||F − Y ||2 = ||U A − Y ||2 .

(a) Montrer que les propositions suivantes sont équivalentes :
– J(A) est un minimum de J sur RN ;
– pour tout B ∈ F, pour tout ε ∈ R, J(A) ≤ J(A + εB);
– pour tout B ∈ F, < U A − Y |U B >= 0.
(b) Vérifier que si M ∈ Mp,q (R), X ∈ Rp , Y ∈ Rq on a

< X|M Y >=<t M X|Y >

où les produits scalaires sont à prendre dans Rq et Rp ...

26
 (c) En déduire que A réalise le minimum de J ssi tU U A −t U Y = 0. Justifier l’exis-
tence et l’unicité de la solution lorsque V est de dimension d + 1.
4. Mise en œuvre sous MAPLE : les données (présentées dans le fichier Moindres-
Carres.mws) sont issues de TIPE de 2004-2005 ;
voir corrigé en 8.0.6

5.2 Algorithme et matrice de Gauss-Newton

1 Pp
Lorsque f (x) = r2 (x), la méthode de Newton appliquée au gradient de f conduit donc
2 k=1 k
à étudier une suite récurrente donnée par Ψ(x) = x − [Dg(x)]−1 g(x) = x − [Hf (x)]−1 g(x) où
−−→ −−→
– g(x) = grad(f )(x) = pk=1 rk (x)grad(rk )(x);
P
– Dg(x) = Hf (x) est la matrice carrée de taille n dont le coefficient de la ligne i colonne j est
p 
∂ 2 rk (x)

X ∂rk (x) ∂rk (x)
+ rk (x) (5.1)
∂xi ∂xj ∂xi ∂xj
k=1
−−→ −−→
Cette matrice est donc la somme de p matrices de rang 1, grad(rk )tgrad(rk ), et des p matrices
rk (x)Hrk (x).  
r1 (x)
• Notons J(x) la matrice jacobienne de x →  ...  , t J(x) est la matrice de Mn,p (R) dont les
 

rp (x)
colonnes sont les p gradients des fonctions rk et le coefficient de la ligne i colonne j de tJ(x) J(x)
est
p
X ∂rk (x) ∂rk (x)
∂xi ∂xj
k=1
On appelle matrice de Gauss-Newton associée problème de minimisation de f la matrice carrée
de taille n, tJ(x) J(x).

Illustration : t J(x) lorsque p = 4, x = 3 :

∂ ∂ ∂ ∂
 
 ∂x1 r1 (x) ∂x1 r2 (x) ∂x1
r3 (x)
∂x1
r4 (x) 
 ∂ ∂ ∂ ∂ 
r1 (x) r2 (x) r3 (x) r4 (x)
 

 ∂x2 ∂x2 ∂x2 ∂x2


 ∂ ∂ ∂ ∂ 
r1 (x) r2 (x) r3 (x) r4 (x)
∂x3 ∂x3 ∂x3 ∂x3
•Comme la k ième colonne de tJ(x) est le gradient de rk (x), nous réécrivons
p
X −−→
g(x) = rk (x) grad(rk )(x) =t J(x)r(x) (5.2)
k=1

Exercice 7 propriétés de la matrice tJ(x) J(x).

1 Pp
On note toujours f (x) = r2 (x) et tJ(x) la matrice de Mn,p (R ) dont les colonnes sont
2 k=1 k
les p gradients des fonctions rk .

27
 1. Montrer que la matrice carrée tJ(x) J(x) est symétrique et positive.
2. Montre que tJ(x) J(x) est définie positive ssi (n ≤ p et J(x) est de rang n).
voir corrigé en 8.0.7

Méthode de Gauss-Newton (Levenberg (40’s) et Marquardt (60’s))

Elle repose sur les idées suivantes :

– On recherche des points stationnaires de f (tels que ||g|| = 0) puisqu’un minimum est nécessairement
atteint en un tel point.
– On adapte la méthode de Newton en remplaçant la matrice hessienne de f par tJ(x) J(x) plus
simple à calculer, à inverser (elle est symétrique positive et on l’espère définie positive). Cela
conduit à poser

xk+1 = xk − λ [tJ(xk ) J(xk )]−1 g(x) = [xk − λ tJ(xk ) J(xk )]−1 tJ(xk )r(xk )

– En choisissant λ réalisant un minimum de la fonction q := λ → f (x − λG−1 gx), nous nous

assurons de plus que (f (xk ))k & .
– Si nous observons la formule (5.1), nous voyons que le second terme y est d’autant plus petit
que les rk y sont proches de 0 (et que les estimations sont bonnes), ou que les dérivées secondes
sont proches de 0 (et que le modèle est proche du linéaire). Il sera peut être possible de prouver
la convergence sous certaines hypothèses.

Algorithme de Gauss-Newton

– On se donne ε > 0, 0 < m1 < m2 < 1 et x1 ;

– On initialise avec
x := x1 ;
−−→
g := grad(f )(x) =t J(x)r(x);
– Tant que ||g|| ≤ ε faire
G :=t J(x)J(x);
q := λ → f (x − λG−1 g);
λ∗ := LineSearchWolfe(q, q1, t0 , m1 , m2 );
(ou tout autre méthode de minimisation en D1)
x := x − λ∗ G−1 g;
−−→
g := grad(f )(x) =t J(x)r(x);
Fin Tant que

Pour implémenter à peu près utilement cet algorithme, il nous faudra une procédure pour le calcul
des fonctions rk , de la fonction f et de son gradient, pour construire les matrices J(x),t J(x) J(x),
la fonction q etc...

28
5.3 Mise en œuvre : MAPLE et formule harmonique
Nous choisissons d’implémenter l’algorithme pour la forme la plus générale de la fonction
p
X
f (x) = rk2 (x)
k=1

Pour cela, avec les notations de la section (5.1), nous choisissons de représenter l’ensemble des
fonctions rk par une fonction du langage dont l’appel suit la syntaxe rk(j, X) où j est un entier,
X la suite des paramètres à estimer, de telle sorte que rk(j, x1 , ..., xn ) = rj (x1 , ..., xn ). Les fonc-
tions à valeurs matricielles associées aux rk , X → J(X) et X →t J(X) J(X) seront construites
formellement par des appels aux procédures MAPLE, J(rk, n, p) et JJ(rk, n, p) comme décrit
dans les tableaux qui suivent. On prendra garde que JJ est une fonction qui retourne une autre
fonction : lorsque r admet la syntaxe précédemment décrite pour rk, l’appel J(r, n, p) retourne
la fonction  
D2 (r) (1, x1 , ..., xn ) . . . Dn+1 (r) (1, x1 , ..., xn )
(x) → 
 .. .. 
. ... . 
D2 (r) (p, x1 , ..., xn ) . . . Dn+1 (r) (p, x1 , ..., xn )
∂
où x est une séquence de n variables et où Di (r)(k, x1 , ..., xn ) = r(k, x1 , ..., xn )
∂xi

Le calcul de J(x) :

J := proc (rk, n, p)
local X, x, i, k, tJ, J;
X := seq(x[i], i = 1 .. n);
seq([seq(diff(rk(k, X),x[i]), i = 1..n)], k = 1..p);
J := Matrix([%]);
unapply(%, seq(X[i], i = 1 .. n))
end proc;

h :=(k, u, v) ->u+vˆ(2*k);
J(h, 2, 3)(a, b);
J(r, 2, 3)(a, b)
 
1 2b
1 4 b3 
1 6 b5
 
D2 (r) (1, a, b) D3 (r) (1, a, b)
D2 (r) (2, a, b) D3 (r) (2, a, b)
D2 (r) (3, a, b) D3 (r) (3, a, b)

29
Le calcul de tJ J(x) : l’appel JJ(r, n, p) retourne la fonction (x) →t J(x) J(x) associée aux
rk (x) = r(k, x).

JJ := proc (rk, n, p)
local X, x, i, k, tJ, J;
X := seq(x[i], i = 1 .. n);
seq([seq(diff(rk(k, X), x[i]), i = 1 .. n)], k = 1 .. p);
J := Matrix([%]);
tJ := Transpose(%);
MatrixMatrixMultiply(tJ, J);
unapply(%, seq(X[i], i = 1 .. n))
end proc;
JJ1:= JJ(v, 2, 3)(x, y):
JJ1[1,2];

D2 (v) (1, x, y) D3 (v) (1, x, y)+D2 (v) (2, x, y) D3 (v) (2, x, y)+D2 (v) (3, x, y) D3 (v) (3, x, y)

Recherche d’un minimum 1D voir justification dans l’exercice 5

GS:=proc(f,a,b,n)
local r,alpha, beta, lambda, mu, k;
r :=(1+sqrt(5))/2;
alpha :=a;
beta :=b;
lambda := evalf(alpha + rˆ2*(beta-alpha));
mu := evalf(alpha + r*(beta-alpha));
for k from 1 to n do
if evalf( f(lambda) - f(mu)) 0
then alpha:= lambda;
else beta := mu;
fi;
lambda := evalf(alpha + rˆ2*(beta-alpha));
mu := evalf(alpha + r*(beta-alpha));
od;
alpha,beta;
end:

30
La méthode sous forme de procédure :
Les arguments :
– rk : fonction (syntaxePrk(j, x1 , ...xn ) = rj (x1 , ..., xn )) ;
– f : fonction (f (x) = pk=1 rk2 (x1 , ...xn ) ; syntaxe f (x1 , ...xn )) ;
– n tel que le nombre d’arguments attendus pour rk est 1 + n; ce sera aussi le nombre de colonnes
de J(rk, n, p);
– p nombre de termes dans la somme des résidus f (x1 , ..., xn ); ce sera aussi le nombre de lignes
de J(rk, n, p);
– eps : flottant donnant un critère d’arrêt ;
– Nmax : nombre maximal d’itérations (si l’on n’a pas atteint auparavant ||g(Xi )|| ≥ eps;
– suite (des arguments : les arguments n˚ 7 à 7 + n sont les termes x01 , ..., x0n de l’initialisation ;
comme le nombre n de ces arguments n’est pas fixe on les appelle avec args[6 + 1], ...args[6 +
n]; on aurait pu choisir de placer là un vecteur mais la gestion des dérivées partielles de rk dans
les procédures J et tJJ aurait été pénible.
−−→
La syntaxe de la valeur retournée : V, Γ (vecteur des paramètres estimé, ||grad(f )(V )||2 )

GaussNewton := proc (rk, f, n, p, eps, Nmax, suite)

local Jr, tJJr, R, g, c, G, iG, V, q, mu, Y;

Jr := J(rk, n, p);
tJJr := tJJ(rk, n, p);

Y := seq(args[6+i], i = 1 .. n);
R := evalf(Vector([seq(rk(k, Y), k = 1 .. p)]));
g := evalf(MatrixVectorMultiply(Transpose(Jr(Y)), R));
c := 0;
while eps<evalf(VectorNorm(g,2,conjugate=false))or Niter <Nmax do
G := evalf(tJJr(Y));
iG := MatrixInverse(G);
V := Vector([Y])- lambda*MatrixVectorMultiply(iG, g);
q := unapply(f(seq(V[i], i = 1 .. 2*N)), lambda);
mu := GS(q, 0, 1, 10);
V := Vector([Y])-mu.MatrixVectorMultiply(iG, g);
Y := seq(V[i], i = 1 .. 2*N);
R := evalf(Vector([seq(rk(k, Y), k = 1 .. P)]));
g := evalf(MatrixVectorMultiply(Transpose(Jr(Y)), R));
c := c+1
end do;

V, evalf(VectorNorm(g, 2, conjugate = false))

end proc;

31
Illustration : amplitudes et phases dans la formule harmonique pour la hauteur de marée

Rappelons brièvement le problème : h0 et les (ωi )1≤i≤N sont connus de même que les relevés
(tk )1≤k≤P et (zk )1≤k≤P (temps et hauteurs mesurés) ; la suite des paramètres à estimer est X =
(a1 , ..., aN , ν1 , ...νN ) on veut donc minimiser
p p
1 X ~ 2 1 X
~ ~ν )2
f (a1 , ..., aN , ν1 , ..., νN ) = φ(A, ~ν , tk ) − zk = rk2 (A,
2 2
k=1 k=1

où les (zk )1≤k≤p sont les hauteurs observées en des instants (tk )1≤k≤p . Nous notons comme dans
notre problème générique a1 , ..., aN , ν1 , ..., νN = x1 , ..., xN , xN +1 , ...x2N . Le nombre de pa-
ramètres est n = 2. le nombre d’observations est p = P.

• Comme toujours, nous allons vérifier nos procédures formellement (débuggage visuel) en expli-
citant les fonctions rk et f puis en testant les procédures J et tJ.
On définit donc en page 33 ces fonctions pour des ωi , νi et des tk , zk formels, ce qui permet de
vérifier et d’illustrer...
• Dans une deuxième étape nous allons tester numériquement nos algorithmes avec des des tk
et zk fictifs pour de petites valeurs de n et p (page 35). Nous verrons déjà apparaı̂tre certains
phénomènes numériques.
• Si tout se passe bien nous passerons à la mise en œuvre sur des données effectives, en l’occur-
rence à partir de 2160 relevés (soit chaque heure pendant 90 jours).

32
Cas particulier des coefficients de marée, représentation formelle des fonctions φ, rk et f
N (on a ici n = 2.N paramètres), P le nombre d’observations, h0 la hauteur moyenne (connue)
0mega = ω1 , ..., ωN suite des pulsations (connues), T suite des instants d’observation et Z suite
des hauteurs observées, sont des variables globales, ici formelles pour illustrer notre propos et
pour permettre de vérifier. X = a1 , ..., aN , ν1 , ..., νN est la suite des variables de φ...

N := 2;
Omega := seq(w[i], i = 1 .. N);

A := seq(a[i], i = 1 .. N);
V := seq(nu[i], i = 1 .. N);
X := A, V;

P := 5;
T := seq(t[i], i = 1 .. P);
Z := seq(z[i], i = 1 .. P);

phi:=(t, X)-> add(args[1+i]*cos(Omega[i]*t+args[N+1+i]),i =1..N);

phi(t, X);

rk :=(k, X)->(phi(T[k], seq(args[j], j = 2..nargs()))-Z[k])ˆ2;

rk(2, X);
rk(4, X);

a1 cos (w1 t2 + ν1 ) + a2 cos (w2 t2 + ν2 ) − z2

a1 cos (w1 t4 + ν1 ) + a2 cos (w2 t4 + ν2 ) − z4

f :=( ) -> (1/2)*add(rk(k, args)ˆ2, k=1..P):

f(X);
1
(a1 cos (w1 t1 + ν1 ) + a2 cos (w2 t1 + ν2 ) − z1 )2
2
1
+ (a1 cos (w1 t2 + ν1 ) + a2 cos (w2 t2 + ν2 ) − z2 )2
2
1
+ (a1 cos (w1 t3 + ν1 ) + a2 cos (w2 t3 + ν2 ) − z3 )2
2
1
+ (a1 cos (w1 t4 + ν1 ) + a2 cos (w2 t4 + ν2 ) − z4 )2
2
1
+ (a1 cos (w1 t5 + ν1 ) + a2 cos (w2 t5 + ν2 ) − z5 )2
2

33
Jr := J(rk, 2*N, P):
Jr(X);
 
cos (w1 t1 + ν1 ) cos (w2 t1 + ν2 ) −a1 sin (w1 t1 + ν1 ) −a2 sin (w2 t1 + ν2 )
cos (w1 t2 + ν1 )
 cos (w2 t2 + ν2 ) −a1 sin (w1 t2 + ν1 ) −a2 sin (w2 t2 + ν2 )
cos (w1 t3 + ν1 )
 cos (w2 t3 + ν2 ) −a1 sin (w1 t3 + ν1 ) −a2 sin (w2 t3 + ν2 )
cos (w1 t4 + ν1 ) cos (w2 t4 + ν2 ) −a1 sin (w1 t4 + ν1 ) −a2 sin (w2 t4 + ν2 )
cos (w1 t5 + ν1 ) cos (w2 t5 + ν2 ) −a1 sin (w1 t5 + ν1 ) −a2 sin (w2 t5 + ν2 )

Test numérique en miniature

• définition des fonctions phi et rk à partir de données numériques simulées
• appel de la procédure (malgré les apparences, le couple (a1 , ν1 ) correspond bien aux valeurs test
car si a1 = −atest
1 , on aν1 = ν1
test + 3π!)

N := 2; P := 5:
Omega := 11, 12:
h0 := 50;
T := seq(k, k = 1 .. P):
phi:=(t,X)->h0+add(args[1+i]*cos(Omega[i]*t+args[N+1+i]),i=1..N);

Xt := 12, 14, 1.6, .7:

Z := NULL;
for k to P do
Z := Z, evalf(phi(T[k], Xt))
end do:

rk :=(k, X)-> phi(T[k], seq(args[j],j = 2..nargs()))-Z[k]:

f :=()->(1/2)*add(rk(k, args)ˆ2, k =1..P):
X0 := 123, 156, 12, 2: eps := 10.ˆ(-6):

X, norme := GaussNewton(rk, f, 2*N, P, eps, 10, X0):

evalf(X[3]-3*Pi);
 
−12.0000000036194461
 14.0000000034084344  −12
 11.0247779606254568  2.420393334 ∗ 10
 

0.699999999879417301
1.59999999862546

34
Test numérique en grandeur nature
• définition des fonctions phi et rk à partir de données numériques observée à raison d’une mesure
toutes les heures pendant 90 jours (p = 24 × 90 = 2160).
• on recherche une évaluation avec n = 8 harmoniques retenues (l’écriture de la fonction f
dépasse 150 pages)
• appel de la procédure

35
5.4 Mise en œuvre : Scilab et formule harmonique

36
6 Méthode quasi-Newton (BFGS)

37
7 Annexe : historique des méthodes de moindres carrés
Pour l’instat, du vrac...
Gauss et Céres
www.scei-concours.fr/tipe/sujet 2009/mathematiques MP 2009.pdf
Cholesky : http ://www.sabix.org/bulletin/b39/vie.html
Toutes les fois que l’on fait une triangulation calculée, il y a avantage à faire également une com-
pensation par le calcul. On est alors amené à écrire un certain nombre d’équations représentant
les relations géométriques entre les divers éléments des figures de la triangulation et comme il
y a généralement plus d’inconnues que d’équations, on lève l’indétermination en écrivant que la
somme des carrés des corrections est minima.
Brezinski http ://smf4.emath.fr/en/Publications/RevueHistoireMath/11/pdf/smf rhm 11 205-238.pdf

38
8 Correction des exercices
Corrigé n˚ 8.0.1 — corrigé de l’exercice 1.
1. C’est la partie facile, on commence par dériver le produit terme à terme, ce qui donne
n
d X d
(A(t) × B(t))i,j = Ai,k (t) × Bk,j (t)
dt dt
k=1
n
X n
X
= A0i,k (t)Bk,j (t) + 0
Ai,k (t)Bk,j (t)
k=1 k=1

d’où la formule de dérivation d’un produit matriciel.

Dans le cas où B est à valeurs dans Rn ,
2. Venons en à l’inversion matricielle J
(a) GLn (R) est l’image réciproque de l’ouvert R∗ par l’application determinant qui est
continue de Mn (R) dans R (les coordonnées sont polynomiales en les mi,j )
(b) J est de classe C ∞ : ce n’est pas encore bien difficile ; on sait que

1 t
A−1 = Com(A)
det(A)

on en déduit que chacune de n2 fonctions composantes de J dans la base canonique

est de classe C ∞ puisque c’est une fonction rationnelle en les mi,j .
(c) Les idées claires.
– pour toute matrice A ∈ GLn (R), dJ (A) est une application linéaire de Mn (R)
dans lui-même.
– Dans la base canonique de Mn (R), la matrice de dJ (A) (matrice jacobienne de
J ) admet donc n2 lignes et n2 colonnes.
∂[A−1 ]k,l
– Chacun des n4 coefficients est de la forme , il s’agit de la dérivée de
ai,j
la composante de la ligne k, colonne l de la fonction A → A−1 par rapport à la
coordonnée ai,j (dans la base canonique).
– Les n2 dérivées partielles de J en un point M dans la base canonique sont les
matrices de Mn (R) définies par

(M + t Ei,j )−1 − M −1 (In + t M Ei,j )−1 − In

 
∂J (A)
= lim = lim M −1
∂mi,j t→0 t t→0 t
3. Venons en au fait :
(a) Pour déterminer dJ (A) qui est donc une application linéaire de Mn (R) dans lui-
même, nous pouvons observer que pour toute courbe paramétrée de classe C 1 , Γ : t ∈
I ∈ R → Γ(t) ∈ GLn (R), on a, par dérivation de fonctions composées,

d(Γ(t))−1 d(J ◦ Γ)
= (t) = dJ (Γ(t)) Γ0 (t)
dt dt

39
(b) On dérive alors la courbe t → Γ(t) × Γ(t)−1 . Il vient avec la formule de dérivation
d’un produit et la formule que précède :

d(Γ(t))−1 × Γ(t) d(J ◦ Γ)(t)) dΓ(t)

= × Γ(t) + J (Γ(t)) ×
dt dt dt
dΓ(t)
= (dJ (Γ(t)) Γ0 (t)) × Γ(t) + +J (Γ(t)) ×
dt
Comme cette expression est nulle (dérivée d’une constante), il vient :

dΓ(t)
dJ (Γ(t)) Γ0 (t) = −Γ(t)−1 × × Γ(t)−1
dt

Cela nous permet de déterminer l’application linéaire dJ (A) pour A ∈ GLn (R). on
considère Γ(t) = A + t H où H est une matrice quelconque, en remplaçant et en
faisant ensuite t = 0, il vient :

dJ (A) H = −A−1 × H × A−1

Remarque : on pourrait aussi déterminer les n2 dérivées partielles de J dans la base

canonique de de Mn (R) en observant que

(M + t Ei,j )−1 − M −1 (In + t M Ei,j )−1 − In

 
∂J (A)
= lim = lim M −1
∂mi,j t→0 t t→0 t

est en fait −M −1 × Ei,j × A−1 ...

40
Corrigé n˚ 8.0.2 — corrigé de l’exercice 2.
1. On dérive donc en adaptant la formule générique de l’exercice 1

dA(t)V (t)
= A0 (t).V (t) + A(t).V 0 (t)
dt
dans laquelle A(t) devient A(X(t)), ... et on obtient par dérivation des fonctions composées
   
d d d
A(X(t))V (X(t)) = A(X(t)) .V (X(t)) + A(X(t)). V (X(t))
dt dt dt
= (dA(X(t)).X 0 (t)).V (X(t)) + A(X(t)).(dV (X(t)).X 0 (t))

Précisons comme l’énoncé nous y invite la nature de chacun des objets que nous avons
introduits :
– dA(X(t)) est une application linéaire de Rn dans Mn (R);
– comme X 0 (t) ∈ Rn , dA(X(t)).X 0 (t) est donc une matrice de Mn (R);
– dA(X(t)).X 0 (t)).V (X(t)) est alors un vecteur de Rn ;
– dV (X(t)) est une application linéaire de Rn dans lui-même ;
– dV (X(t)).X 0 (t) ∈ Rn de même que A(X(t)).(dV (X(t)).X 0 (t)).
2. Considérons alors la courbe X : t → x0 + t.h, notons A.V l’application x → A(x).V (x).
Il vient, en dérivant la fonction composée (A.V ) ◦ X et en rapprochant de ce qui précède :

d
A(X(t))V (X(t)) = d(A.V )(X(t)).X 0 (t)
dt
d
A(X(t))V (X(t)) = (dA(X(t)).X 0 (t)).V (X(t)) + A(X(t)).(dV (X(t)).X 0 (t))
dt
ce qui, lorsque t = 0, donne

d(A.V )(x0 )h = dA(x0 ).h.V (x0 ) + A(x0 ).V (x0 ).h

41
Corrigé n˚ 8.0.3 — corrigé de l’exercice 3

1. La fonction f est de classe C 2 , φ : x → det df (x) est donc de classe C 1 et F est définie là
où ce déterminant ne s’annule pas. Or l’image réciproque de R∗ par φ est un ouvert contenu
dans Ω.
Par ailleurs, en utilisant en partie le résultat de l’exercice 2, nous avons

F (x) = x − (J ◦ df )(x).f (x)

dF (x) = In − {d(J ◦ df ))(x).f (x) + (J ◦ df )(x).df (x)}
= In − {d(J ◦ df ))(x).f (x) + In }

Sans aller plus loin dans les calculs, avec f (c) = 0 on a F (c) = c et dF (c) = 0 ∈ L(Rn ).
2. Comme dF (c) = 0 et comme dF est continue, il existe δ1 > 0 tel que sur B1 = B̄(c, δ1 )
|||dF (x)||| ≤ 1.
(a) Pour x = c + h ∈ B1 , nous avons
Z 1
d
F (c + h) − F (c) = F (c + t h) dt
0 dt
Z 1
= dF (c + t h) h dt
0
Z 1
||F (c + h) − F (c)|| = |||dF (c + t h)||| ||h|| dt dt
0

Ainsi, ||F (c + h) − c|| ≤ ||h|| ≤ δ1 , et cette boule est bien stable par F.
(b) On reprend la démonstration vue en dimension 1, on écrit

xn+1 − c = xn − c − [Df (xn )]−1 .f (xn ) (8.1)

−1
= −[Df (xn )] (f (xn ) + Df (xn )(c − xn )) (8.2)

Ecrivons la formule de Taylor à l’ordre 2 au voisinage d’un point x ∈ B,

1
f (c) =h→0 f (x) + Df (x).(c − x) + D2 f (x)(c − x, c − x) + ε(c − x) ||(c − x)||2
2

Elle permet d’écrire sur la boule B1 privée de c,

−||c − x||2
   
2 c−x c−x
f (x) + Df (x).(c − x) = D f (x) , + ~ε(c − x)
h→0 2 ||c − x|| ||c − x||
   
2 c−x c−x
** L’expression D f (x) , + ~ε(c − x) est bornée sur B1 privée
||c − x|| ||c − x||
de c. En effet,
– (x, h, h) → D2 f (x)(h, h) est continue ; l’expression donnée page 11 le montre
clairement.

42
– Elle est donc bornée sur B1 × Sn−1 × Sn−1 qui est un compact de Rn × Rn × Rn .
– Enfin, x → ~ε(c−x) a pour limite 0 en c. Elle est donc bornée sur une boule B̄(0, δ2 ).
On note alors B = B̄(0, min(δ1 , δ2 ) = B̄(0, δ). Cette boule est encore stable par F
(puisque |||DF (x)||| y est encore majorée par 1) et, en reprenant la formule (8.1) pour
une suite non stationnaire,

 ||c − xn ||2
   
−1
 2 c − xn c − xn
||xn+1 − c|| = |||Df (xn )] ||| D f (x)
 , + ~ε(c − xn )
||c − xn || ||c − xn || 2

ce qui conduit à la majoration attendue avec

1
K= max |||Df (x)]−1 || × max ||D2 f (x)(h, h) + ~ε(h)||
2 x∈B x∈B,||h||=1

43
Corrigé n˚ 8.0.4 — corrigé de l’exercice 4.
Soit f : I → R une fonction définie sur un intervalle I de R.
1. Nous allons procéder en deux étapes avant de conclure.
• Étape 1 : On suppose tout d’abord que x1 < x2 = x3 < x4 . La relation du théorème
devient :

f (x2 ) − f (x1 ) f (x4 ) − f (x2 )

≤ . (8.3)
x2 − x1 x4 − x2
Comme les dénominateurs sont strictement positifs, on multiplie par chacun d’eux ce qui
nous donne la relation équivalente :

(f (x2 ) − f (x1 ))(x4 − x2 ) ≤ (f (x4 ) − f (x2 ))(x2 − x1 )

elle-même équivalente à
   
x4 − x2 x2 − x1
f (x2 ) ≤ f (x1 ) + f (x4 ).
x4 − x1 x4 − x1
Il s’agit là de la définition même de la convexité de f puisque
   
x4 − x2 x2 − x1
x2 = x1 + x4 .
x4 − x1 x4 − x1

Cela montre que f est convexe ssi (8.3).

• Étape 2 :
Si (x1 , x2 , x2 , x4 ) ∈ I 4 est tel que x1 < x2 ≤ x3 < x4 , on a en appliquant deux fois le
résultat (8.3)

f (x2 ) − f (x1 ) f (x3 ) − f (x2 ) f (x4 ) − f (x3 )

≤ ≤ . (8.4)
x2 − x1 x3 − x2 x4 − x3
Comme (8.3) est un cas particulier de (8.4), on a aussi (8.3)⇔ (8.4).

• Ainsi f convexe ⇔ (8.3)⇔ (8.4).

2. On suppose f est strictement convexe et continue sur I = [a, b]. Elle y admet un minimum
puisqu’elle est continue sur un segment (compact).
Montrons que ce minimum est atteint en un point au plus en raisonnant par l’absurde : du
fait de la stricte convexité, si ce minimum est atteint en deux points distincts, x1 et x2 avec

f (x1 ) = f (x2 ) = min f (x) = µ,

x∈[a,b]

on aura en particulier,
 
1 1 1 1
f x1 + x2 < f (x1 ) + f (x2 ) = µ.
2 2 2 2
Cela contredit la définition du minimum.

44
3. On se donne t ∈]0, 1/2[ et on partage I = [a, b] en trois intervalles de la façon suivante :

a+b

x = − t (b − a)
 1

 2
I = [a, x1 ] ∪ [x1 , x2 ] ∪ [x2 , b] avec
x2 = a + b + t (b − a)



2
(a) On note p0 , p2 , p2 les pentes de f sur ces trois intervalles.
– Montrons que p0 ≤ p1 ≤ p2 ≤ 0 =⇒ µ ∈ [x2 , b];
On suppose donc que

f (x1 ) − f (a) f (x2 ) − f (x1 ) f (b) − f (x2 )

p0 = < p1 = < p2 = < 0. (8.5)
x1 − a x2 − x1 b − x2
f (b) − f (x2 )
Comme p2 = ≤ 0, d’après le résultat de la première question, sur
b − x2
[a, x2 ] le taux de variation de f est strictement négatif et f décroı̂t strictement sur
cet intervalle. Donc, si x ∈ [a, x2 ],

f (x) > f (x2 ) ≥ min f (t)

t∈[a,b]

le minimum n’est donc pas atteint sur [a, x2 ], il l’est nécessairement sur [x2 , b].

– Montrons que p0 ≤ p1 < 0 =⇒ µ ∈ [x1 , b].

Toujours d’après le résultat de la première question, comme p1 < 0, sur [a, x1 ]
le taux de variation de f est strictement négatif, f décroı̂t sur x ∈ [a, x1 ] et si
x ∈ [a, x1 ]
f (x) > f (x1 ) ≥ min f (t)
t∈[a,b]

de cela on déduit que le minimum est atteint sur [x1 , b] puisqu’il ne peut l’être sur
[a, x1 ].
– Montrons que p0 < 0 < p1 < p2 =⇒ µ ∈ [a, x2 ].
Comme p2 > 0, sur [x2 , b] le taux de variation de f est strictement positif et f est
strictement croissante. Donc, pour x ∈ [x2 , b], on a :

min f (t) ≤ f (x2 ) < f (x).

t∈[a,b]

Le minimum ne saurait être atteint sur [x2 , b] il est donc atteint sur [a, x2 ].
– on montre de la même façon que 0 ≤ p0 ≤ p2 ≤ p2 =⇒ µ ∈ [a, x1 ]...
(b) Voir pour cela le chapitre 6 de [1].

45
Corrigé n˚ 8.0.5 — corrigé de l’exercice 5
1. I = [a, b] et f & strict. sur ] − ∞, t0 ] ∩ I et f % strict. sur [t0 , +∞[∩I.
√
5−1
(a) Observons que si r est la solution positive de x2
+ x = 1, r = ∈]0, 1[,
√ 2
3− 5
r2 = 1 − r = ∈]0, 1[, ce qui assure que an ≤ λn ≤ µn ≤ bn lorsque
2
λn = an + r2 `n et µn = an + r`n , avec `n = bn − an .

Nous allons montrer par récurrence sur n que t0 ∈ [an , bn ].

C’est évident lorsque n = 0. Supposons que t0 ∈ [an , bn ]. Deux cas se présentant
alors :
• Soit f (λn ) ≥ f (µn ).
Cela exclut que t0 ∈ [an , λn [, sinon, f strictement croissante sur [t0 , bn ] on aurait
f (t0 ) ≤ f (λn ) < f (µn ) ce qui est incompatible avec l’hypothèses.
Comme on pose dans ce cas an+1 = λn et bn+1 = bn la relation t0 ∈ [an+1 , bn+1 ] est
établie.
• Si f (λn ) ≤ f (µn ) on ne peut avoir t0 ∈ [µn , bn ] pour une raison analogue et comme
an+1 = an et bn+1 = µn la relation t0 ∈ [an+1 , bn+1 ] est établie.
Convergence :
on a `n+1 ≤ max{r, r2 } = r`n donc `n ≤ rn `0 et les suites (an )n croissante, et (bn )n
décroissante sont adjacentes et encadrent t0 ..
(b) GS:=proc(f,a,b,n)
global r;
local alpha,beta,lambda,mu, k;

alpha :=a;
beta :=b;

lambda := evalf(alpha + rˆ2*(beta-alpha));

mu := evalf(alpha + r*(beta-alpha));

for k from 1 to n do
if evalf( f(lambda) - f(mu)) 0
then alpha:= lambda;
else beta := mu;
fi;
lambda := evalf(alpha + rˆ2*(beta-alpha));
mu := evalf(alpha + r*(beta-alpha));
od;

alpha,beta;
end:

46
 g:=proc(x)
if x<=1 then 3-x
elif x <=2 then x+1
else 3;
fi;
end;

plot(g,0..3,0..3);
fsolve(3*rˆn=10ˆ(-5));
GS(g,0,3,27);
g := proc (x) if x <= 1 then 3−x elif x <= 2 then x+1 else 3 end if end proc

26.20787166
0.9999971839, 1.000004014

2. Minimum de fonctions de plusieurs variables

F est définie sur un ouvert U convexe de R2 , à valeurs réelles, de classe C 1 .
(a)

dφ(λ) ∂F (x + λu) d(x1 + λu1 ) ∂F (x + λu) d(x2 + λu2 )

= + =< grad(F )(x+λu)|u >
dλ ∂x1 dλ ∂x2 dλ
(règle de la chaı̂ne)
Ainsi, φ0 (0) =< grad(F )(x)|u >= ||grad(F )(x)||.||u||. cos(θ), expression dans la-
grad(F )(x)
quelle seule θ est variable : φ0 (0) est minimale lorsque θ = π et u = − .
||grad(F )(x)||
(b) On définit la suite avec
• si grad(F )(xn ) 6= 0,
−grad(F )(xn )
– un =
||grad(f )(xn )||
– xn+1 = xn + λ∗n un où λ∗n réalise le minimum de F (xn + λun ).

47
 • si grad(F )(xn ) = 0, xn est le dernier terme de la suite qui est donc finie.

F (xn+1 ) = F (xn + λ∗n un ) ≤ F (xn + 0un ) = F (xn ), la suite (F (xn ))n est donc
décroissante et minorée (puisque F est minorée) et elle est convergente.

(c) Exemple F (x, y) = ln 1 + x2 + y 2 + xy .

f:=(x,y)->ln(1+xˆ2+yˆ2+x*y);
plot3d(f(x,y),x=-1..1,y=-1..1);

[diff(f(x,y),x),diff(f(x,y),y)];
g:=unapply(%,(x,y));

2x + y 2y + x
(x, y) 7→ [ , ]
1 + x2 + y 2 + xy 1 + x2 + y 2 + xy
n:=12;
X:=[1.,3];

for k from 1 to n do
u := normalize(g(X[1],X[2]));
phi := lambda - f( X[1] + lambda*u[1], X[2] + lambda*u[2]);
m := GS(phi,-5,5,15);
lambdaStar := (m[1]+m[2])/2;
X := evalm(X+lambdaStar*u);
f(X[1],X[2]);
od;
n := 12
X := [1., 3]
u := [0.5812381936, 0.8137334712]

48
φ := λ− > f (X[1] + λu[1], X[2] + λu[2])
m := −2.921626370, −2.914294996
lambdaStar := −2.917960683
X := [−0.696030196, 0.625557725]
0.3649024098
..
.
X := [−0.001860518352, 0.0008530456106]
0.260299661210−5

49
Corrigé n˚ 8.0.6 corrigé de l’exercice 6 Approximation discrète au sens des moindres carrés
1. Généralités
(a) Si f est combinaison linéaire des d + 1 fonctions (ui )0≤i≤d , il existe (ai )i ∈ Kd+1 tel
que f (x) = di=0 ai ui (x) et alors
P

   Pd   
f (x1 ) i=0 ai ui (x1 ) d ui (x1 )
F =  ...  =  ..  X  .. 
= ai  .  .
  
.
Pd i=0
f (xN ) i=0 ai ui (xN )
ui (xN )
 
ui (x1 )
Ainsi le vecteur F est combinaison linéaire des d + 1 vecteurs  ...  . Nous
 

ui (xN )
noterons V le sous-espace de K N engendré par ces d+1 vecteurs. On peut exprimer
matriciellement :
 
u0 (x1 ) u1 (x1 ) . . . ud (x1 )  
 u0 (x2 ) u1 (x2 ) . . . ud (x2 )  a0
 
  a1 
F = UA =  .
  
. .
. .
.   .. 
 . . .   .
 

  ad
u0 (xN ) u1 (xN ) . . . ud (xN )

avec U ∈ MN,d+1 et A ∈ Kd+1 . V apparaı̂t alors comme étant l’image de U dans

K N . Cette image a pour dimension le rang de U qui est au plus d + 1.
(b) Lorsque la fonction f est polynomiale de degré ≤ d, elle est combinaison linéaire
des fonctions monômes ui (x) = xi et la matrice U est la matrice de Vandermonde
(tronquée) :

1 x1 . . . xd1
 
 
u0 (x1 ) u1 (x1 ) . . . ud (x1 ) d
 u0 (x2 ) u1 (x2 ) . . . ud (x2 )  1 x2 . . . x2 

   
   
U = . = .
   
.. ..  
 .. . .   .
. .
. .
. 
  . . . 
   
 
u0 (xN ) u1 (xN ) . . . ud (xN ) d
1 xN . . . x N

2. Considérons une norme || || dans KN . D’après les préliminaires topologiques, pour tout
Y ∈ K N , il existe un vecteur F0 ∈ V (non nécessairement unique) de V sev de KN tel que

||Y − F0 || = inf ||Y − F ||.

F ∈V

3. On pose J(A) = ||F − Y ||2 = ||U A − Y ||2 .

(a) Numérotons :
i. J(A) est un minimum de J sur RN ;

50
 ii. pour tout B ∈ F, pour tout ε ∈ R, J(A) ≤ J(A + εB);
iii. pour tout B ∈ F, < U A − Y |U B >= 0.
(1) ⇒ (2) est immédiat.
(2) ⇒ (3) : la fonction de la variable réelle ε → J(A + εB) = ||U (A + εB) − Y ||22
s’exprime :

J(A + εB) = ||U (A + εB) − Y ||22 = ||U B||2 ε2 + 2ε < U B|U A − Y > +J(A).

Si pour tout ε, J(A) ≤ J(A + εB), on a aussi pour tout ε > 0,

||U B||2 ε + 2 < U B|U A − Y >≥ 0.

Cela impose (par exemple en passant à la limite) < U B|U A − Y >≥ 0.

De la même façon, si ε < 0,

||U B||2 ε + 2 < U B|U A − Y >≤ 0,

ce qui impose < U B|U A − Y >≤ 0.

(3) ⇒ (1) : on suppose que pour tout B ∈ KN , < U B|U A − Y >= 0. On a alors
pour X ∈ Kd+1 quelconque

J(X) = J(A+(X−A)) = ||U (X−A)||2 +2 < U (X−A)|U A−Y > +J(A) ≥ J(A).

(b) Soient M ∈ Mp,q (R), X ∈ Rp , Y ∈ Rq ; on a M Y ∈ Rp et

p q p q
! !
X X X X
< X|M Y >= xi mi,k yk = mi,k xi yk .
i=1 k=1 i=1 k=1

q p q p
! !
X X X X
t
< M X|Y >= yk (t m)k,i xi = mi,k xi yk
k=1 i=1 k=1 i=1

où les deux expressions sont égales.

(c) La relation < U A − Y |U B >= 0 qui devient <t U (U A − Y )|B >= 0. Comme B
est quelconque, F = U A réalise le minimum de ||F − Y || ssi tU U A −t U Y = 0.
Une telle équation admet une solution A et une seule lorsque tU U est inversible. Or
la matrice carrée tU U a le même noyau que U. Elle est inversible dans Md+1 (R) dès
lors que rg(U ) = dimV = d + 1.

Vérifions que ker(U ) = ker(t U U ).

U ∈ MN,d+1 (R),t U U ∈ Md+1 (R) considérons X ∈ Rd + 1.
-X ∈ ker(U ) implique que pour tout Y, < U X|U Y >=<t U U X|Y >= 0. Ainsi
t U U X = 0.

- Réciproquement, si t U U X = 0, <t U U X|X >=< U X|U X >= ||U X||2 = 0 et

X ∈ ker(U ).
Remarque : l’existence est assurée quelque soit le rang de U ou la dimension de V
pour des raisons topologiques.

51
Corrigé n˚ 8.0.7 corrigé de l’exercice 7
1 Pp
On note toujours f (x) = r2 (x) et tJ(x) la matrice de Mn,p (R ) dont les colonnes sont
2 k=1 k
les p gradients des fonctions rk .
1. Montrer que la matrice carrée tJ(x) J(x) est symétrique et positive.
2. Montre que tJ(x) J(x) est définie positive ssi (n ≤ p et J(x) est de rang n).

52
Références
[1] T HIERRY AUDIBERT, A MAR O USSALAH
Informatique, Programmation et Calcul Scientifique
Ellipses, 2013
[2] M. B ERGONIOUX .
Optimisation et contrôle des systèmes linéaires
Dunod, 2001
[3] J.F R ÉD ÉRIC B ONNANS , J.C HARLES G ILBERT, C LAUDE L EMAR ÉCHAL , C LAUDIA A. S A -
GASTZ ÁBAL .
Numerical Optimisation
Springer, 2006
[4] C LAUDE B REZINSKI , M ICHELA R EDIVO -Z AGLIA .
Méthodes numériques itératives
Ellipse, 2006
[5] J.P. D EMAILLY .
Analyse numérique et équations différentielles
PUG, 1996
[6] R ICHARD A. H OLMGREN .
A first course in dynamical systems
Springer, 2000
[7] SHOM .
http ://www.shom.fr/fr page/fr act oceano/maree/ondes.htm
la page du SHOM, Service hydrographique de la Marine
[8] B ERNARD S IMON .
La marée 1. La Marée océanique cotière
Institut océanographique Collection : SYNTHESES (2007)
[9] JAN A. S NYMAN .
Practical Mathematical Optimisation
Springer, 2005

53
Index
attractif moindres carrés
point fixe, 6 problème, 24

bassin d’attraction, 5, 22 Newton

quasi (méthodes), 16
Cholesky norme
moindres carrés, 38 subordonnée, 14
convergence
géométrique, 5 ordre de convergence, 5
linéaire, 5 méthode de Newton, 7
quadratique, 5 méthode du point fixe, 6
vitesse de, 5
point fixe
différentielle, 10, 11 attractif, 6
de l’inversion matricielle, 12 super attractif, 2
seconde, 11
régression
formule linéaire
de Taylor, 10, 11 une variable, 26
règle
Gauss-Newton d’or (rech. minimum), 18
matrice de, 27
système
hessienne, 10
dynamique, 5
jacobienne, 10
Taylor
ligne de descente formule de, 10, 11
règle d’or, 30
Wolfe
règle de Wolfe, 19
procédé (rech. minimum), 19
méthode
quasi Newton, 16
méthode de Newton
plusieurs variables, 14
une variable, 2
matrice
de Gauss-Newton, 27
hessienne, 10
jacobienne, 10
minimum
fonction convexe, 17
ligne de descente
règle d’or, 30
règle de Wolfe, 19
règle d’or, 18

54