Module : Optimisation Année 2020-2021

Solution de la série de TD 4

Exercice 1 La méthode de plus forte pente est utilisée pour résoudre le problème

min f (x) = 2x21 − 2x1 x2 + x22 + 2x1 − 2x2

et une séquence {xk } est générée.
1. En supposant que   T
1
x2k+1 = 0 1− k
5
Montrer que   T
1
x2k+3 = 0 1 − k+1
5
2. Trouver un minimiseur de f (x) en utilisant les résultats de (1).

Solution. Soient
min f (x) = 2x21 − 2x1 x2 + x22 + 2x1 − 2x2
et une séquence {xk } est générée.

1. La méthode de la plus forte pente

gkT gk
xk+1 = xk − gk
gkT Hk gk

où g est le gradient et H est le hessien.

   
4x1 − 2x2 + 2 4 −2
g(x) = et H(x) =
2x2 − 2x1 − 2 −2 2

Dans cet exercice changeons l'indice k par 2k + 1, ainsi la formule devient :

T
g2k+1 g2k+1
x2k+2 = x2k+1 − T
g2k+1
g2k+1 H2k+1 g2k+1
En supposant que

  T
1
x2k+1 = 0 1− k
5

1
nous avons  
2  
 5k  4 −2
g2k+1 (x) =  et H2k+1 (x) = 
  

 2  −2 2
− k
5
Après calcul, on a :

−2
 
 5k+1 
x2k+2 =
 

 3 
1−
5k+1
Après le même raisonnement on trouve :

 
0
x2k+3 = 
 
1

1−
5k+1
C.Q.F.D.

%Déclaration de la fonction

function [y]=fct(x1,x2)
y = 2 ∗ x1 − 2 ∗ x1 ∗ x2 + x22 + 2 ∗ x1 − 2 ∗ x2 ;
2

end

%Calcul du gradient

function [G]=grad(x)
x1 = x(1) ;
x2 = x(2) ;
G = [4 ∗ x1 − 2 ∗ x2 + 2; −2 ∗ x1 + 2 ∗ x2 − 2] ;

%Calcul du hessien

function [H]=hess(x)
H = [4 -2 ;-2 2] ;

%Calcul de x2k+3

function [ak,bk] = exo1source(fname,gname,hname,x0)

xk = x0 ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
ak = xk-((gk'*gk)/(gk'*Hk*gk))*gk ;
gkk = feval(gname,ak) ;
Hkk = feval(hname,ak) ;
bk=ak-((gkk'*gkk)/(gkk'*Hkk*gkk))*gkk ;
disp('x2k+2=')

2
 ak
disp('x2k+3=')
bk

Pour avoir les résultats souhaités, exécuter dans la command window

syms k

0
[ak, bk] = exo1source('fct','grad','hess', 0 1 − 1/5k )

2. En faisant tendre k vers ∞ on obtient x? = (0 1)T

Exercice 2 Le problème
min f (x) = x21 + 2x22 + 4x1 + 4x2
est résolu en utilisant la méthode de la plus forte pente avec comme point initial x1 = (0 0)T .
1. Par récurrence, montrer que
"  k #T
2 1
xk+1 = k −2 − −1
3 3

2. Déduire le minimiseur de f (x).

Solution.

1. Par récurrence, montrer que

"  k #T
2 1
xk+1 = k −2 − −1
3 3

g1T g1
Tout d'abord, calculons x2 à partir de la formule x2 = x1 − g1 . Pour cela calculons
g1T H1 g1
le gradient et le hessien de cette fonction.

   
2x1 + 4 2 0
g(x) = et H(x) =
4x2 + 4 0 4
 
4  
4


4 4
 −3 
     
0 4 4
x2 =  −      = 



0 2 0 4 4
 4 
−


4 4    
3
0 4 4
Vérions que cette formule est vraie pour k = 1.
"  1 #T  T
2 1 4 4
x2 = 1 − 2 − −1 = − −
3 3 3 3
Maintenant, supposons que cette formule est vraie pour k =n−1 c'est-à-dire

"  n−1 #T
2 1
xn = −2 − −1
3n−1 3

3
 et montrons qu'elle resteras vraie pour k=n c'est-à-dire

  n T
2 1
xn+1 = n −2 − −1
3 3

Pour cela calculons xn+1 à partir de la formule de la méthode de la plus forte pente
donnée dans l'exercice 1.

gnT gn
xn+1 = xn − gn
gnT Hn gn

donc
 4 
  n−1   3n−1
4 1

4 −
 
 2 
3n−1 3
 n−1 
1
 4 
−2
 
3n−1 4 − 3n−1
  3  
xn+1 =   −
   
 1
n−1 
 4    n−1
 1


− −1 4 −
 
  n−1  2 0  3n−1
3 4 1
 3
4 − 
 
  n−1 
3n−1 3 0 4
 1 
4 −
3
Ce qui donne

 2   4 
−2
 3n−1  1 3n−1 
xn+1 =   − 
   
 3
n−1  n−1 
 1 1 
− −1 4 −
3 3
Ainsi  2 
−2
 3n 
xn+1 = 
 
 n 
 1 
− −1
3
C.Q.F.D.

2. En faisant tendre k vers ∞, on trouve x? = (−2 − 1)T .

Exercice 3 Soit le problème

min f (x) = x21 + x22 − 0.2x1 x2 − 2.2x1 + 2.2x2 + 2.2

1. Trouver un point vériant les conditions nécessaires du premier ordre pour un minimi-
seur.
2. Montrer que ce minimiseur est global.
3. Quel est le taux de convergence de la méthode de plus forte pente pour ce problème ?
4. En partant du point x0 = [0 0]T , combien d'itérations de plus forte pente au plus sont
elles nécessaires pour réduire la valeur de la fonction à 10−10 .

4
Solution.
Soit
min f (x) = x21 + x22 − 0.2x1 x2 − 2.2x1 + 2.2x2 + 2.2
1.  
2x1 − 0.2x2 − 2.2
∇f (x) = g(x) =  
2x2 − 0.2x1 + 2.2
Conditions du premier ordre g(x? ) = 0. Ainsi le point vériant ces conditions c'est
x? = (1, −1).
2.

3. Le taux de convergence de la méthode de la plus forte pente.

f (xk+1 ) − f (x? ) ≤ β [f (xk ) − f (x? )]

où β est le taux de convergence déni par

 2
1−r
β=
1+r
avec
la plus petite valeur de Hk
r=
la plus grande valeur de Hk

Vérions tout d'abord que H(x? ) est bien déni positif. Donc

 
2 −0.2
H(x? ) =  
−0.2 2
Ici, D = 3.96 > 0 donc H est déni positif. Maintenant, calculons les valeurs propres de
H.

2 − λ −0.2
P (λ) = |H − λI| = = (2 − λ)2 − (0.2)2 = (1.8 − λ)(2.2 − λ)
−0.2 2 − λ
P (λ) = 0 ⇐⇒ λ1 = 1.8 et λ2 = 2.2
ainsi,
1.8
r= = 0.8182
2.2
par suite
 2
1 − 0.8182
β= = 0.0099
1 + 0.8182
4. Nous avons f (1, −1) = 0 et f (0, 0) = 2.2. ainsi

 2  4  2k
1−r 1−r 1−r
f (xk ) ≤ f (xk−1 ) ≤ f (xk−2 ) ≤ · · · ≤ f (0, 0)
1+r 1+r 1+r | {z }
2.2

et on résout l'équation
 2k
−10 1−r
10 = 2.2
1+r

5
 ce qui donne  
1−r
−10 ln(10) = ln(2.2) + 2k ln
1+r
Donc le nombre d'itérations est : k = 2.
%Déclaration de la fonction

function [y]=fct(x1,x2)
2 2
y = x1 + x2 − 0.2 ∗ x1 ∗ x2 − 2.2 ∗ x1 + 2.2 ∗ x2 + 2.2 ;
end

%Calcul du gradient

function [G]=grad(x)
x1 = x(1) ;
x2 = x(2) ;
G = [2 ∗ x1 − 0.2 ∗ x2 − 2.2; 2 ∗ x2 − 0.2 ∗ x1 + 2.2] ;

%Calcul du hessien

function [H]=hess(x)
H = [2 -0.2 ;-0.2 2] ;

%Algorithme de la méthode de la plus forte pente sans recherche linéaire avec hessien

function [xs,fs,k] = steep_desc(fname,gname,hname,x0,epsi)

disp(' ')
disp('Program steep_desc.m')
k = 1;
xk = x0 ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
ak = (gk'*gk)/(gk'*Hk*gk) ;
adk = -ak*gk ;
er = norm(adk) ;
while er >= epsi,
xk = xk + adk ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
ak = (gk'*gk)/(gk'*Hk*gk) ;
adk = -ak*gk ;
er = norm(adk) ;
k = k + 1;
end
format long
disp('Solution point :')

6
 xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

Résultats
[xs,fs,k] = steep_desc('fct','grad','hess',[0 0]',1e-10)
Program steep_desc.m
Solution point :
xs =(1 -1)
Objective function at the solution point :
fs =0
Number of iterations performed :
k = 2

Exercice 4 Soit
min f (x) = 5x21 − 9x1 x2 + 4.075x22 + x1
1. Résoudre ce problème en utilisant la méthode de la plus forte pente avec x0 = [1 1]T et
 = 3.10−6 .
2. Faire une analyse de convergence du problème pour expliquer pourquoi la méthode de la
plus forte pente requiert ici un grand nombre d'itérations pour atteindre la solution.

Solution.

1. %Déclaration de la fonction

function [y]=fct(x1,x2)
y = 5 ∗ x1 − 9 ∗ x1 ∗ x2 + 4.075 ∗ x22 + x1 ;
2

end

%Calcul du gradient

function [G]=grad(x)
x1 = x(1) ;
x2 = x(2) ;
G = [10 ∗ x1 − 9 ∗ x2 + 1; 8.15 ∗ x2 − 9 ∗ x1] ;

%Calcul du hessien

function [H]=hess(x)
H = [10 -9 ;-9 8.15] ;

En utilisant l'algorithme de la méthode de la plus forte pente sans recherche linéaire

avec hessien, on trouve comme résultats :

[xs,fs,k] = steep_desc('fct','grad','hess',[1 1]',3*1e-6)

Program steep_desc.m

7
 Solution point :
xs =(-16.299620676374516 -17.999579247889098)
Objective function at the solution point :
fs =-8.149999995571662
Number of iterations performed :
k =1347

2. Analyse de convergence.

Calculons le taux de convergence déni par

 2
1−r
β=
1+r
avec
la plus petite valeur de Hk
r=
la plus grande valeur de Hk

Vérions tout d'abord que H(x? ) est bien déni positif. Donc

 
10 −9
H(x? ) =  
−9 8.15
Ici, D = 0.5 > 0 donc H est bien déni positif. Maintenant, calculons les valeurs propres
de H.

10 − λ −9
P (λ) = |H − λI| = = (10 − λ) (8.15 − λ) − 81
−9 8.15 − λ
P (λ) = 0 ⇐⇒ λ1 = 18.12 et λ2 = 0.02
ainsi,
18.12
r= = 0.0015
0.02
par suite
 2
1 − 0.0015
β= = 0.9939
1 + 0.0015
Exercice 5 Résoudre le problème
min f (x) = 100(x2 − x21 )2 + (1 − x1 )2
en utilisant
1. l'algorithme de la plus forte pente sans recherche linéaire avec hessien.
2. l'algorithme de la plus forte pente sans recherche linéaire sans hessien.
3. l'algorithme de la plus forte pente avec recherche linéaire et vérier la solution en utilisant
les conditions susantes de second ordre.
4. l'algorithme de Newton.
5. l'algorithme de Gauss-Newton.
6. la commande prédénie de MATLAB fminunc
dans les cas suivants
x0 = [4 4]T , x0 = [4 − 4]T , x0 = [−4 4]T et x0 = [−4 − 4]T avec  = 10−6 .
Comparer les résultats obtenus.

8
 Solution.

%Déclaration de la fonction

function f = fonction(x)
x1 = x(1);
x2 = x(2);
f = 100 ∗ (x2 − x12 )2 + (1 − x1)2 ;

%Déclaration du gradient

function g = gradient(x)
x1 = x(1) ;
x2 = x(2) ;
g1 = −400 ∗ x1 ∗ (x2 − x12 ) − 2 ∗ (1 − x1) ;
g2 = 200 ∗ (x2 − x12 ) ;
g = [g1 g2]0 ;

%Déclaration du hessien

function H = hessien(x)
x1 = x(1);
x2 = x(2);
h11 = 1200 ∗ x12 − 400 ∗ x2 + 2 ;
h12 = −400 ∗ x1 ;
h22 = 200 ;
H = [h11 h12; h12 h22] ;

%Déclaration du jacobien

function J = jacobien(x)
x1 = x(1);
x2 = x(2);
J = [−20 ∗ x1 10; −1 0] ;

1. Méthode 1 : Algorithme de la plus forte pente sans recherche linéaire avec

hessien

function [xs,fs,k] = steep_desc1(fname,gname,hname,x0,epsi)

disp(' ')
disp('Program steep_desc1.m')
k = 1;
xk = x0 ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
ak = (gk'*gk)/(gk'*Hk*gk) ;

9
 adk = -ak*gk ;
er = norm(adk) ;
while er >= epsi,
xk = xk + adk ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
ak = (gk'*gk)/(gk'*Hk*gk) ;
adk = -ak*gk ;
er = norm(adk) ;
k = k + 1;
end
format long
disp('Solution point :')
xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

2. Méthode 2 : Algorithme de la plus forte pente sans recherche linéaire sans

hessien

function [xs,fs,k] = steep_desc2(fname,gname,x0,epsi)

disp(' ')
disp('Program steep_desc2.m')
k = 1;
xk = x0 ;
fk = feval(fname,xk) ;
gk = feval(gname,xk) ;
ah = 1 ;
fh = feval(fname,xk-ah*gk) ;
gk2 = gk'*gk ;
2
ak = (gk2*ah )/(2*(fh-fk+ah*gk2)) ;
adk = -ak*gk ;
er = norm(adk) ;
while er >= epsi,
xk = xk + adk ;
fk = feval(fname,xk) ;
gk = feval(gname,xk) ;
ah = ak ;
fh = feval(fname,xk-ah*gk) ;
gk2 = gk'*gk ;
2
ak = (gk2*ah )/(2*(fh-fk+ah*gk2)) ;
adk = -ak*gk ;
er = norm(adk) ;
k = k + 1;
end
format long

10
 disp('Solution point :')
xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

3. Méthode 3 : Algorithme de la plus forte pente avec recherche linéaire

function [xs,fs,k] = steep_desc3(fname,gname,x0,epsi)

disp(' ')
disp('Program steep_desc3.m')
k = 1;
xk = x0 ;
gk = feval(gname,xk) ;
dk = -gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = norm(adk) ;
while er >= epsi,
xk = xk + adk ;
gk = feval(gname,xk) ;
dk = -gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = norm(adk) ;
k = k + 1;
end
format short
disp('Solution point :')
xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

Algorithme de la recherche linéaire.

function z = inex_lsearch(xk,s,F,G,p1,p2)
k = 0;
m = 0;
tau = 0.1 ;
chi = 0.75 ;
rho = 0.1 ;
sigma = 0.1 ;
mhat = 400 ;
epsilon = 1e-10 ;

11
xk = xk( :) ;
s = s( :) ;
parameterstring = ;
% evaluate given parameters :
if nargin > 4,
if isstr(p1),
eval([p1 ' ;']) ;
else
parameterstring = ',p1' ;
end
end
if nargin > 5,
if isstr(p2),
eval([p2 ' ;']) ;
else
parameterstring = ',p2' ;
end
end
% compute f0 and g0
eval(['f0 = ' F '(xk' parameterstring ') ;']) ;
eval(['gk = ' G '(xk' parameterstring ') ;']) ;
m = m+2 ;
deltaf0 = f0 ;
% step 2 Initialize line search
dk = s ;
aL = 0 ;
aU = 1e99 ;
fL = f0 ;
dfL = gk'*dk ;
if abs(dfL) > epsilon,
a0 = -2*deltaf0/dfL ;
else
a0 = 1 ;
end
if ((a0 <= 1e-9)|(a0 > 1)),
a0 = 1 ;
end
% step 3
while 1,
deltak = a0*dk ;
eval(['f0 = ' F '(xk+deltak' parameterstring ') ;']) ;
m = m + 1;
% step 4
if ((f0 > (fL + rho*(a0 - aL)*dfL)) & (abs(fL - f0) > epsilon) & (m < mhat))
if (a0 < aU)
aU = a0 ;
end
% step 5 compute a0hat using equation 7.65
2
a0hat = aL + ((a0 - aL) *dfL)/(2*(fL - f0 + (a0 - aL)*dfL)) ;

12
 a0Lhat = aL + tau*(aU - aL) ;
if (a0hat < a0Lhat)
a0hat = a0Lhat ;
end
a0Uhat = aU - tau*(aU - aL) ;
if (a0hat > a0Uhat)
a0hat = a0Uhat ;
end
a0 = a0hat ;
else
eval(['gtemp =' G '(xk+a0*dk' parameterstring ') ;']) ;
df0 = gtemp'*dk ;
m = m + 1;
% step 6
if (((df0 < sigma*dfL) & (abs(fL - f0) > epsilon) & (m < mhat) & (dfL ∼= df0)))
deltaa0 = (a0 - aL)*df0/(dfL - df0) ;
if (deltaa0 <= 0)
a0hat = 2*a0 ;
else a0hat = a0 + deltaa0 ;
end
a0Uhat = a0 + chi*(aU - a0) ;
if (a0hat > a0Uhat)
a0hat = a0Uhat ;
end
aL = a0 ;
a0 = a0hat ;
fL = f0 ;
dfL = df0 ;
else
break ;
end
end
end% while 1
if a0 < 1e-5,
z = 1e-5 ;
else
z = a0 ;
end

4. Méthode 4 : Algorithme de Newton

function [xs,fs,k] = newton(fname,gname,hname,x0,dt,epsi)

disp(' ')
disp('Program newton.m')
k = 1;
n = length(x0) ;
In = eye(n,n) ;
xk = x0 ;
gk = feval(gname,xk) ;

13
 Hk = feval(hname,xk) ;
[V, D] = eig(Hk) ;
di = diag(D) ;
dmin = min(di) ;
if dmin > 0,
Hki = V*diag(1./di)*V' ;
else
bt = dt - dmin ;
Hki = V*diag((1+bt)./(di+bt))*V' ;
end
dk = -Hki*gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = norm(adk) ;
while er >= epsi,
xk = xk + adk ;
gk = feval(gname,xk) ;
Hk = feval(hname,xk) ;
[V, D] = eig(Hk) ;
di = diag(D) ;
dmin = min(di) ;
if dmin > 0,
Hki = V*diag(1./di)*V' ;
else
bt = dt - dmin ;
Hki = V*diag((1+bt)./(di+bt))*V' ;
end
dk = -Hki*gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = norm(adk) ;
k = k + 1;
end
format long
disp('Solution point :')
xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

5. Méthode 5 : Algorithme de Gauss-Newton

function [xs,fs,k] = gauss_newton(fname,gname,jname,x0,epsi)

disp(' ')
disp('Program gauss_newton.m')
x = x0( :) ;
n = length(x) ;

14
 In = eye(n) ;
k = 1;
xk = x0 ;
F_k = feval(fname,xk) ;
gk = feval(gname,xk) ;
Jk = feval(jname,xk) ;
Hk = 2*Jk'*Jk + 1e-12*In ;
dk = -inv(Hk)*gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = norm(adk) ;
while er > epsi,
xk = xk + adk ;
F_k1 = feval(fname,xk) ;
gk = feval(gname,xk) ;
Jk = feval(jname,xk) ;
Hk = 2*Jk'*Jk + 1e-12*In ;
dk = -inv(Hk)*gk ;
ak = inex_lsearch(xk,dk,fname,gname) ;
adk = ak*dk ;
er = abs(F_k1 - F_k) ;
k = k + 1;
F_k = F_k1 ;
end
format long
disp('Solution point :')
xs = xk + adk
disp('Objective function at the solution point :')
fs = feval(fname,xs)
format short
disp('Number of iterations performed :')
k

6. La commande fminunc
fminunc solutionne des problèmes d'optimisation non linéaires et multi-variables et sans
restrictions. Cette fonction permet de changer entre algorithmes diérents par exemple
de méthode de Interior Reective Newton (si on connait les dérivées) ou le Méthode
BFGS ( Broyden, Fletcher, Goldfarb et Shanno) en cas contraire. L'algorithme BFGS
approxime la matrice Hessienne (méthode quasi-Newton).

[x,fval,outout,exitag] = fminunc('f ',x0)

Le minimum de cette fonction est : x? = (1, 1).

15
 Tableau récapitulatif.

Conditions Initiales [4, 4]T [−4, 4]T [4, −4]T [−4, −4]T

Nbre d'itérations de la méthode 1 5687 4358 269 10199

Nbre d'itérations de la méthode 2 23672 11740 9085 9103

Nbre d'itérations de la méthode 3 15148 8948 10733 11089

Nbre d'itérations de la méthode 4 16 20 13 21

Nbre d'itérations de la méthode 5 4 6 4 5

fminunc 32 47 41 42

On remarque dans cette exemple que la méthode de Gauss-Newton est une méthode op-
timale. Comme les autres méthodes numériques étudiées dans cet exemple, elle est également
basée sur des calculs itératifs. Bien que les propriétés théoriques de convergence soient très
bonnes (convergence quadratique au voisinage de la solution).
Remarque. Toutes ces méthodes numériques ont pour conséquences pratiques de dénir
un "point d'initialisation".
Les calculs vont commencer en ce point dans l'espace de recherche. Le choix de ce point
d'initialisation ne doit pas être négligé. Parfois un mauvais choix est synonyme d'échec de la
mise en oeuvre de la méthode.

16