Zangré Ibrahim

Optimisation non linéaire sans contraintes

Rapport du TP2

Exercice 1 : Minimisation d’une fonctionnelle non quadratique

On définit sur R2 la fonctionnelle de Rosenbrock par :

f (x, y) = (x − 1)2 + 10(x2 − y)2 (1)

1. Etude théorique
(a) Points critiques de f :
f est C ∞ et on a :
 
2(x − 1) + 40x(x2 − y)
∇f (x, y) =  
2
−20(x − y)
et :
 
120x2 − 40y + 2 −40x
Hess [f ] (x, y) =  
−40x 20
On a : 
 (x − 1) + 20x(x2 − y) = 0
∇f (x, y) = 0 ⇐⇒
x2 − y = 0



 x−1=0
⇐⇒
x2 − y = 0



 x=1
⇐⇒

y=1  
1
Le seul point critique de f est donc X ∗ :=
1

1
 (b)
On a f (X ∗ ) = 0 et ∀X ∈ R2 , f (X) ≤ f (X ∗ ) = 0 ; donc X ∗ est un minimum global de f ; comme c’est
l’unique point critique de f , alors c’est l’unique minimum global.

 
∗ 82 −40
(c) Soit H la matrice hessienne de f calculé au point X ; H = ; son polynôme caracté-
−40 20
√ √
51− 2561 51+ 2561
ristique est P = λ2 −102λ+40, dont les racines sont : λ1 = 2 (= λmin ) et λ2 = 2 (= λmax ) ;
le conditionnement de H est donné par :
λmax
κ(A) = ≈ 258, 096
λmin

(d)Devéloppment
 à l’ordre 2 de f en X ∗ :
u1
soit u := , alors :
u2
   
1
f (X ∗ + u) = f (X ∗ ) + ∇f (X ∗ ), u + Hu, u + kuk2 ǫ(u) = 41u21 + 10u22 − 40u1 u2 + kuk2 ǫ(u)
| {z } | {z } 2
0 ~0

2. Etude numérique

(a) Représentation de la fonctionnelle sur le pavé [−10, 10] × [−10, 10] :

Fonction Matlab :

function []=plotGraphe(h);

% h : pas de la discrétisation spaciale

x=-10:h:10;
y=-10:h:10;
[x,y]=meshgrid(x,y);
z=(x-1).^2+10.*(x.^2-y).^2;

surf(x,y,z);

Ce qui donne :

2
 4
x 10

14

12

10

2
10
0
−10 0
−5
0
5
10 −10

Figure 1 – Représentation sur le pavé [−10, 10] × [−10, 10]

(b) Tracer de lignes de niveau de la fonctionnelle

Fonction Matlab :

function []=plot_ligne_de_nivau(a,delta,N);

% a :borne du pavé
% delta :pas de discrétisation
% N :nombre de lignes de niveau
% à afficher

x=-a:delta:a;
y=-a:delta:a;
[x,y]=meshgrid(x,y);
z=(x-1).^2+10.*(x.^2-y).^2;

[c,h]=contour(x,y,z,N);
colormap hsv

Ce qui donne :

3
 10

−2

−4

−6

−8

−10
−10 −5 0 5 10

Figure 2 – Lignes de niveau

La “platitude” de la surface au voisinage du point X ∗ ralentit la convergence des algorithmes : à partir

d’un certain moment “on avance plus” !
3. Comparaison de deux méthodes de minimisation
(a) Méthode du gradient à pas optimal

i. Minimisation de la fonction α → f (X − α.∇f (X)) avec la fonction Matlab f minsearch :

function [alphaOptimal,fval]=pasOptimal(alpha0,X);

[dfx,dfy]=gradientF(X);

rosenbrok=@(alpha)((X(1)-alpha.*dfx)-1).^2+10.*((X(1)-alpha.*dfx).^2-(X(2)-alpha.*dfy)).^2;

[alphaOptimal,fval] = fminsearch(rosenbrok,alpha0);

ii, iii. Gradient à pas optimal

On utilisera un double critère d’arrêt de l’algorithme :
-un critère portant sur la précision de la méthode, par exemple portant sur le module du gradient
-et un critère portant sur le nombre d’itérations pour éviter que l’algorithme soit infini.

4
function [X_etoile]=AlgoPasOtptimal(X0,N,tol);

[dfx,dfy]=gradientF(X0);
X(:,1)=X0;
test=norm([dfx,dfy],2);
k=1;

while test>tol && k<=N

[dfx,dfy]=gradientF(X(:,k));
[alphaOptimal,fval]=pasOptimal(1,X(:,k));
X(:,k+1)=X(:,k)-alphaOptimal.*[dfx;dfy];
%X(k,:)=X(k-1,:)-0.5/(51+sqrt(2561)).*[dx2 dy2];

test=norm([dfx,dfy],2);

k=k+1;
end

X_etoile=X(:,k-1);

figure
hold on
plot_ligne_de_niveau(10,0.5,50)
plot(X(1,:),X(2,:),’g’);
plot(1,1,’*r’)
hold off

   
0 0, 9999
En partant du point X0 = ,avec N = 103 et tol = 10−3 , on obtient : X ∗ = : on représente
1 0, 9998
les lignes de niveau et la trajectoires des x(k) :

5
 2

1.5

0.5

−0.5

−1

−1.5

−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

Figure 3 – Lignes de niveau et trajectoire des x(k)

(b) Méthode du gradient à pas fixe

Fonction Matlab : pasF ixe

function [X_etoile]=AlgoPasFixe(X0,beta,N,tol);

[dfx,dfy]=gradientF(X0);
X(:,1)=X0;
test=norm([dfx,dfy],2);
k=1;

while test>tol && k<=N

[dfx,dfy]=gradientF(X(:,k));
X(:,k+1)=X(:,k)-beta.*[dfx;dfy];

test=norm([dfx,dfy],2);

k=k+1;
end

X_etoile=X(:,k-1);

Xop=AlgoPasOtptimal([0;1],1e+2,1e-3);
figure

6
 hold on
plotGraphe2(2,0.1,50)
plot(X(1,:),X(2,:),’r’,’LineWidth’,2)
legend(’pas fixe’);
plot(Xop(1,:),Xop(2,:),’color’,[0 1 1])
legend(’pas optimal’);
plot(1,1,’*b’)
hold off

 
3 −3 ∗ 0, 9989
Avec N = 10 et tol = 10 , on obtient : X = .
0, 9977
(c) Comparaison des méthodes

i. Trajectoires des x(k) générés par les gradients à pas optimal et à pas fixe :

1.5

0.5

−0.5

−1

−1.5

−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

Figure 4 – Pas optimal (vert) et pas fixe (rouge)

ii. Les lignes de niveau sont rares au voisinage de X ∗ ; cela est dû à la “platitude” de la fonc-
tionnelle au voisinage de ce point ; la valeur du pas donné par la recherche 1D optimale est à l’origine
du zizag des x(k) générés par le gradient à pas optimal (la courbe est plus lisse pour des petits pas : voir
gradient à pas fixe → fig4).

(k)
iii. On modifie les programmes prédédents pour qu’ils revoient laliste des points
 x obtenus
(k) ∗
par les deux algorithmes ; puis on calcule et on trace la courbe donnant ln kx − X k en fonction de

7
k;

Fonction Matlab :

function []=erreur();

Xop=AlgoPasOtptimal([0;1],1e+3,1e-3);
Xfix=AlgoPasFixe([0;1],0.01,1e+3,1e-3);
N_op=size(Xop,2);
N_fix=size(Xfix,2);
X_etoile=[1;1];

for i=1:N_op
erreur_op(i)=log(norm(Xop(:,i)-X_etoile,2));
end

for i=1:N_fix
erreur_fix(i)=log(norm(Xfix(:,i)-X_etoile,2));
end
hold on
plot((1:N_op),erreur_op,’g’)
plot((1:N_fix),erreur_fix,’r’)
hold off

soit :

1
pas optimal
pas fixe
0

−1

−2

−3

−4

−5

−6

−7
0 200 400 600 800 1000 1200

 
(k) ∗
Figure 5 – Fonction erreur ln kx −X k

Ces courbes ont une allure linéaire ; l’erreur kx(k) − X ∗ k peut être donnée par e−ck , avec c > 0. On voit

8
bien que la gradient conjugué converge plus vite que le gradient à pas fixe (pente plus grande en module).

Exercice 2 : Résolution de systèmes linéaires à l’aide de préconditioneurs

On souhaite résoudre le système d’inconnue x : An x = bn , avec :

 
2 −1 0 ... 0
 ..   
−1 2 −1 . 
  1

An =  0 .. .. .. 
∈ Mn (R), et bn = . . . ∈ Rn
 . . . 0 
 .  1
 .. −1 2 −1
0 ... 0 −1 2

1. Première méthode : le gradient conjugué.

Résoudre le système précédent revient à minimisé la fonctionnelle : Jn (x) = 21 (An x, x) − (bn , x)

Algorithme du gradient conjugué


x(0) donné, r(0) = An x(0) − b, d(0) = −r(0)

 (n) (n)
= x(n) + ρn d(n) , ρn = − (A(rn d(n)
,d )
 (n+1)
 x ,d(n) )



 r(n+1) = An x(n+1) − b


 (n+1) 2
d(n+1) = −r(n+1) + βn d(n) , βn = krkr(n) k2k

Focntion Matlab : Gradient conjugué

function [x_etoile,k]=gradientConjugue(n,x0,Niter,tol);

%Niter : nombre maximal d’itérations

%tol : test d’arrêt sur le module du gradient
%k : donne le nombre d’itérations nécéssaires pour
% atteindre la solution
A=An(n);
x(:,1)=x0;
r(:,1)=dJn(n,x(:,1));
d(:,1)=-r(:,1);
i=1;
k=0;
while test>tol && i<=Niter-1

k=k+1;

rho=-((d(:,i))’*r(:,i))/((d(:,i))’*(A*d(:,i)));
x(:,i+1)=x(:,i)+rho*d(:,i);

9
 r(:,i+1)=dJn(n,x(:,i+1));
beta=(norm(r(:,i+1),2)/norm(r(:,i),2))^2;
d(:,i+1)=-r(:,i+1)+beta*d(:,i);

test=norm(d(:,i),2);

i=i+1;

end

x_etoile=x(:,i-1);

2. La méthode du gradient conjugué préconditionné.

Si M est une matrice inversible, le système An x = bn est identique à : M −1 An x = M −1 bn ; on choisit la

matrice M pour que le conditionnement de M −1 An soit meilleur que celui de An , et que M −1 soit facile
à calculer.

(a) Principe de la factorisation de Cholesky

Soit A une matrice symétrique définie positive ; alors il existe une unique matrice triangulaire à élé-
ments diagonaux positifs R telle que :
t
A = RR
C’est un cas particulier de la factorisation LU : soit à résoudre le système Ax = b ; ce système est équi-
valent à :
 t
 Ry = b

Rx = y
(b) La factorisation peut “remplir” partiellement la matrice R ; donc on considère RI la factorisée
incomplète de Cholesky de An , et on pose M = t RI × RI ; dans ce cas la matrice M est facile à inverser ;
comme M est “proche” de A, M −1 A est proche de I, et son conditionnement plus proche de 1, donc
meilleur que celui de A .

(c) Gradient conjugué préconditionné :

function [Xtt,k]=GradientConjuguePreconditionne(n,X0,Niter,tol);

%Niter : nombre maximal d’itérations

%tol : test d’arrêt sur le module du gradient
%k : donne le nombre d’itérations nécéssaires pour
% atteindre la solution

An=MatrixCreator(n);
RI=cholinc(sparse(An),’0’);

10
M=RI’*RI;
Minv=inv(M);

x(:,1)=X0;
r(:,1)=Minv*dJn(n,x(:,1));
d(:,1)=-r(:,1);
test=norm(d(:,1),2);
i=1;
k=0;
while test>tol && i<=Niter-1

k=k+1;

rho=-((d(:,i))’*r(:,i))/((d(:,i))’*(Minv*An*d(:,i)));
x(:,i+1)=x(:,i)+rho*d(:,i);
r(:,i+1)=Minv*dJn(n,x(:,i+1));
beta=(norm(r(:,i+1),2)/norm(r(:,i),2))^2;
d(:,i+1)=-r(:,i+1)+beta*d(:,i);

test=norm(d(:,i),2);

i=i+1;

end

X_etoile=x(:,i-1);
b=ones(n,1);
X_theoric=An\b;
Xtt=[X_etoile,X_theoric];

On améliore ainsi la vitesse de convergence ; la solution est atteinte en deux itérations !

3. Comparaison des algorithmes.

On complète le programme pour qu’il renvoie aussi le résidu logarithmique en fonction du nombre d’ité-
rations ;

Le résidu logarithmique pour la méthode du gradient conjugué préconditionné est à −∞ en quelques

itérations : méthode nettement plus rapide que le gradient conjugué sans préconditionnement.

11
 5 grad conj
grad conj précond

−5

−10

−15

−20

−25

−30

−35
0 200 400 600 800 1000

Figure 6 – Résidu logarithmique

12