Résolution De Système Linéaire - Méthode Itérative Gradient Conjugué - JULIA


Ce petit “notebook” a pour but de montrer quelques propriétés sur l’algorithme du gradient conjugué préconditionné. Le langage utilisé est Julia, l’interface utilisée -le notebook- est Ipython avec le package IJulia. Les graphiques sont réalisés à l’aide du package PyPlot.

Les différents tests (benchmarks) seront menés sur un système d’équation provenant de la discrétisation différence finie du problème modèle :

$$\begin{cases}
-\Delta u = f & \text{dans } \Omega \
u=0& \text{sur } \partial \Omega
\end{cases}$$

Le 11/04/2016

Thierry Clopeau (clopeau@univ-lyon1.fr)


1. Assemblage du système

voir https://forgeim.univ-lyon1.fr/blog/2016/03/15/SysLin_Methode_Directe/

1
using PyPlot # bibliothèque graphique
1
2
3
4
5
6
function Laplace1d(n)
h=1./(n+1)
d1=fill(-1/h^2,n-1);
d2=fill(2/h^2,n);
spdiagm((d1,d2,d1),(-1,0,1));
end
Laplace1d (generic function with 1 method)
1
2
3
function Laplace2d(n,m)
kron(speye(m),Laplace1d(n))+kron(Laplace1d(m),speye(n))
end
Laplace2d (generic function with 1 method)
1
2
3
function Laplace3d(n,m,p)
kron(speye(p),Laplace2d(n,m))+kron(Laplace1d(p),speye(n*m))
end
Laplace3d (generic function with 1 method)

2. Algorithme du Gradient Conjugué préconditionné

L’algorithme :

  • Precond est une fonction passée en paramètre, on peut passer la fonction identité avec x->x en paramètre i.e. pcg(A,b,x->x,x0,1e-12,1000).
  • Le test d’arrêt est calculé sous la forme $\frac{|r|}{|b|} \leq \varepsilon$ i.e. le résidu normalisé par le second membre.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function pcg(A,b,Precond,x0,Eps,itermax)
# Renormalisation et test
bnrm=norm(b);
if bnrm==0
return zeros(b)
end
#initialisation
x=x0;
r=b-A*x;
z=Precond(r);
d=copy(z);
res=[norm(r)/bnrm];
tmp=dot(r,z);
iter=1;

while (res[iter]>Eps) & (iter<itermax)
Ad=A*d;
alpha=tmp/dot(d,Ad);
x=x+alpha*d;
r=r-alpha*Ad;
z=Precond(r);
tmp1=dot(r,z);
beta=tmp1/tmp;
d=z+beta*d;
tmp=tmp1;
iter=iter+1;
push!(res,norm(r)/bnrm);
end
return x,res
end
pcg (generic function with 1 method)

Préconditionnement - préambule

En premier lieu on définit 2 fonctions de résolution (efficace) de sytèmes triangulaires (inf et sup). Ces 2 fonctions remplacent aventageusement celles utilisées par la fonction \ de JULIA. Les matrices creuses sont en format CSC.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function triusolve(A,b)
if !istriu(A) throw(DimensionMismatch("A is'nt upper triangular matrix")) end
n = length(b)
if isa(b, Vector)
nrowB = n
ncolB = 1
else
nrowB, ncolB = size(b)
end
ncol = sqrt(prod(size(A)))
if nrowB != ncol throw(DimensionMismatch("A is $(ncol)X$(ncol) and b has length $(n)")) end

x=zeros(n);
for i=n:-1:1
x[i]=(b[i]-x[i])/A[i,i];
for j=A.colptr[i]:A.colptr[i+1]-2;
x[A.rowval[j]]+=A.nzval[j]*x[i]
end
end
return x
end
triusolve (generic function with 1 method)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function trilsolve(A,b)
if !istril(A) throw(DimensionMismatch("A is'nt lower triangular matrix")) end
n = length(b)
if isa(b, Vector)
nrowB = n
ncolB = 1
else
nrowB, ncolB = size(b)
end
ncol = sqrt(prod(size(A)))
if nrowB != ncol throw(DimensionMismatch("A is $(ncol)X$(ncol) and b has length $(n)")) end

x=zeros(n);
for i=1:n
x[i]=(b[i]-x[i])/A[i,i];
for j=A.colptr[i]+1:A.colptr[i+1]-1;
x[A.rowval[j]]+=A.nzval[j]*x[i]
end
end
return x
end
trilsolve (generic function with 1 method)

Préconditionnement - Diagonal

On utilise une variable global $D$ afin de stocker la diagonale de la matrice.

Le calcul se fait en 2 étapes, dans un premier temps on assemble le préconditionneur (ici $D$) puis à chaque itération on procède à la résolution du système donné par le préconditionneur. Le choix en terme de coût et d’inverser 1 fois la diagonale de la matrice et ensuite de multiplier à chaque itération (la division coûte plus chère que la multiplication).

Attention : dans notre cas la digonale est de la forme $\lambda I_d$ et donc le préconditionnement diagonal n’améliore pas du tout le conditionnement du système !

1
2
3
4
5
6
7
8
9
10
function init_PrecDiag(A)
global D
D=1 ./ diag(A);
return "PrecDiag initialised"
end

function PrecDiag(r)
global D
z= D .* r;
end
PrecDiag (generic function with 1 method)
1
2
3
4
function PrecDiag(r)
global D
z= D .* r;
end
PrecDiag (generic function with 1 method)

Préconditionnement - SSOR

On utilise les variable global $U$, $L$, $D$ et $w$ afin de stocker la décomposition de la matrice de préconditionnement. L’assemblage n’est fait qu’une seule fois :
$$ M_\omega= \frac{\omega}{2-\omega}\left(\frac{D}{\omega}+E\right)\,D^{-1}\,\left(\frac{D}{\omega}+E^T\right) \quad E\text{ partie strictement supérieure de }A$$

d’où

$$ \left(M_\omega\right)^{-1}= \left(\frac{D}{\omega}+E^T\right)^{-1} \left(\frac{2-\omega}{\omega} D \right)\left(\frac{D}{\omega}+E\right)^{-1}$$

Le calcul se fait en 2 étapes, dans un premier temps on assemble le préconditionneur :

  • $L \leftarrow \left(\frac{D}{\omega}+E\right)$
  • $D \leftarrow \frac{2-\omega}{\omega} D$
  • $U\leftarrow\left(\frac{D}{\omega}+E^T\right)$.

A chaque itération on procède à la résolution du système à l’aide des 2 fonctions précédentes triusolve et trilsolve fonction de résolution de système triangulaires u-pper et l-ower.

Pour une formule de $w_{opt}$ voir http://userpages.umbc.edu/~gobbert/papers/YangGobbert2007SOR.pdf

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function init_SSOR(A,w_in)
global U,L,D,w
w=copy(w_in)
D=diag(A)/w;
U=triu(A,1)+spdiagm(D);
L=tril(A,-1)+spdiagm(D);
D=D/(2-w);
return "SSOR initialised"
end

function PrecSSOR(r)
global U,L,D,w
r=trilsolve(L,r);
r=D.*r
r=triusolve(U,r);
end
PrecSSOR (generic function with 1 method)

Evolution du résidu

Regardons pour une matrice de taille respectable (Laplace1d, 2d et 3d) l’évolution du résidu (normalisé) par $b$, regardons l’effet des préconditionneurs sur la décroissance du résidu nous utiliserons la valeur optimal pour SSOr : $w_{opt}=\frac{2}{1+\sin(\pi h)}$.

Résidu - Laplace1d

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
nn=1000000;
A=Laplace1d(nn);
b=ones(nn);
h=1/(nn+1);
wopt=2/(1+sin(pi*h));

## Sans préconditionnement
(x,res1)=pcg(A,b,x->x,zeros(b),1e-12,400);
## préconditionnement diagonal
init_PrecDiag(A)
(x,res2)=pcg(A,b,PrecDiag,zeros(b),1e-12,400);
## Préconditionnement SSOR
init_SSOR(A,wopt)
(x,res3)=pcg(A,b,PrecSSOR,zeros(b),1e-18,400);
f=figure()
xlabel("Itération k")
ylabel("Norme du résidu")
title("Evolution du résidu - Laplace 1d de taille $nn")
semilogy(1:length(res1),res1,color="red")
semilogy(1:length(res2),res2,color="blue")
semilogy(1:length(res3),res3,color="green")
legend(["Sans préconditionneur" ; "Précond. Diagonal" ;"Précond. SSOR"],loc=4)
grid("on")

png

Constation

  • Les itérations sans préconditionnement semblent ne pas converger ou très très doucement…
  • Le préconditionnement diagonal est sans effet (par rapport au non préconditionnement). Cela est normal car la diagonale est de la forme une constante fois l’identité, on ne change pas le conditionnement d’une matrice en a multipliant par une constante !
  • Le péconditionnement SSOR est efficace, seulement 18 itérations pour converger vers un résidu en $10^{-16}$ pour une matrice de taille $10^6$.

Résidu Laplace2d

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
nn=1000;
A=Laplace2d(nn,nn);
b=ones(nn^2);
h=1/(nn+1);
wopt=2/(1+sin(pi*h));

## Sans préconditionnement
(x,res1)=pcg(A,b,x->x,zeros(b),1e-12,400);
## préconditionnement diagonal
init_PrecDiag(A)
(x,res2)=pcg(A,b,PrecDiag,zeros(b),1e-16,400);
## Préconditionnement SSOR
init_SSOR(A,wopt)
(x,res3)=pcg(A,b,PrecSSOR,zeros(b),1e-18,400);

figure();
xlabel("Itération k")
ylabel("Norme du résidu")
title("Evolution du résidu - Laplace 2d de taille $(nn*nn)")
semilogy(1:length(res1),res1,color="red")
semilogy(1:length(res2),res2,color="blue")
semilogy(1:length(res3),res3,color="green")
legend(["Sans préconditionneur" ; "Précond. Diagonal" ;"Précond. SSOR"],loc=5)
grid("on")

png

Constation

  • La convergence sans préconditionnement est plus marquée qu’en 1d mais pas extra-ordinaire !
  • Le préconditionnement diagonal est toujours sans effet…
  • Le péconditionnement SSOR est peu efficace, avec seulement plus de 250 itérations pour converger vers un résidu en $10^{-16}$ pour une matrice de taille $10^6$.

Résidu Laplace3d

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
nn=100;
A=Laplace3d(nn,nn,nn);
b=ones(nn^3);
h=1/(nn+1);
wopt=2/(1+sin(pi*h));

## Sans préconditionnement
(x,res1)=pcg(A,b,x->x,zeros(b),1e-12,400);
## préconditionnement diagonal
init_PrecDiag(A)
(x,res2)=pcg(A,b,PrecDiag,zeros(b),1e-12,400);
## Préconditionnement SSOR
init_SSOR(A,wopt)
(x,res3)=pcg(A,b,PrecSSOR,zeros(b),1e-18,400);

figure()
xlabel("Itération k")
ylabel("Norme du résidu")
title("Evolution du résidu - Laplace 3d de taille $(nn^3)")
semilogy(1:length(res1),res1,color="red")
semilogy(1:length(res2),res2,color="blue")
semilogy(1:length(res3),res3,color="green")
legend(["Sans préconditionneur" ; "Précond. Diagonal" ;"Précond. SSOR"],loc=4)
grid("on")

png

Constation

  • La convergence sans préconditionnement est encore plus marquée qu’en 1d et 2d rendant l’algoritme compétitif !
  • Le préconditionnement diagonal est toujours sans effet…
  • Le péconditionnement SSOR est efficace, avec seulement 75 itérations pour converger vers un résidu en $10^{-16}$ pour une matrice de taille $10^6$.

3. Performance du Gradient Conjugué préconditionné (SSOR)

Nous venons de voir que la bibliothèque UMFPACK possède des limites vites atteintes (particulierement en 3d). Le gradient conjugué avec un préconditionnement SSOR est une alternative, elle permet au moins dans le cas 3d d’atteindre des taille de matrices plus importantes, le côut mémoire étant fortement réduit. Le préconditionnement SSOR réduit fortement le nombre d’itérations du Gradient conjugué.

Mais le gradient conjugué préconditionné SSOR est-il vraiment plus performant que la bibliothèque UMFPACK ?

Il faut donc comparer les temps de résolution, en prenant également en compte le temps de l’assemblage du préconditionneur !

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# pour enlever le temps de compilation de la fonction \
Laplace1d(3)\ones(3);
# début
npt=12
println("1d... ");
N1d=round(Int,logspace(2,7,npt));
T1d=zeros(length(N1d));
T1d1=zeros(length(N1d));
for i=1:length(N1d)
A=Laplace1d(N1d[i]); b=ones(N1d[i]);
wopt=2/(1+sin(pi/(N1d[i]+1)));
println(" n=$(N1d[i])")
# opérateur \
tic();A\b;T1d[i]=toq();
# GC preconditionné
tic();
init_SSOR(A,wopt)
(x,rr)=pcg(A,b,PrecSSOR,zeros(b),1e-16*N1d[i],400);
T1d1[i]=toq();
end;
println(" completed!");

println("2d... ");
N2d=round(Int,logspace(1,3.5001,npt));
T2d=zeros(length(N2d));
T2d1=zeros(length(N2d));
for i=1:length(N2d)
A=Laplace2d(N2d[i],N2d[i]);b=ones(N2d[i].^2)
wopt=2/(1+sin(pi/(N1d[i]+1)));
println(" n=$(N2d[i])")
# opérateur \
tic();A\b;T2d[i]=toq();
# GC preconditionné
tic();
init_SSOR(A,wopt)
(x,rr)=pcg(A,b,PrecSSOR,zeros(b),1e-16*N2d[i]^2,400);
T2d1[i]=toq();
end;
println(" completed!");

println("3d...");
N3d1=round(Int,logspace(1,2.3334,npt));
N3d=N3d1[1:8]
T3d=zeros(length(N3d));
T3d1=zeros(length(N3d1));
for i=1:length(N3d1)
A=Laplace3d(N3d1[i],N3d1[i],N3d1[i]);b=ones(N3d1[i].^3)
wopt=2/(1+sin(pi/(N1d[i]+1)));
println(" n=$(N3d1[i])")
# opérateur \ - limité en taille maximale
if (i<=8)
tic();A\b;T3d[i]=toq();
end
# GC preconditionné
tic();
init_SSOR(A,wopt)
(x,rr)=pcg(A,b,PrecSSOR,zeros(b),1e-16*N3d1[i].^3,400);
T3d1[i]=toq();
end;
println(" completed!");
1d... 
  n=100
  n=285
  n=811
  n=2310
  n=6579
  n=18738
  n=53367
  n=151991
  n=432876
  n=1232847
  n=3511192
  n=10000000
           completed!
2d... 
  n=10
  n=17
  n=28
  n=48
  n=81
  n=137
  n=231
  n=390
  n=658
  n=1111
  n=1874
  n=3163
          completed!
3d...
  n=10
  n=13
  n=17
  n=23
  n=31
  n=40
  n=53
  n=71
  n=93
  n=123
  n=163
  n=215
          completed!
1
2
3
4
5
6
7
8
9
10
11
12
figure()
xlabel("taille de la matrice")
ylabel("temps en s")
title("Temps de résolution UMFPACK - GC SSOR")
loglog(N1d,T1d,color="red")
loglog(N1d,T1d1,color="red",linestyle="--")
loglog(N2d.^2,T2d,color="blue")
loglog(N2d.^2,T2d1,color="blue",linestyle="--")
loglog(N3d.^3,T3d,color="green")
loglog(N3d1.^3,T3d1,color="green",linestyle="--")
legend(["1d umfpack"; "1d pgc" ; "2d umfpack" ;"2d pgc" ;"3d umfpack" ;"3d pgc"],loc=2)
grid("on")

png

Analyse et conclusion

Regardons ces dernières courbes :

  • On voit que UMFPACK côute cher à l’allumage, il est peu compétitf sur de petites tailles de matrice…

Pour des matrices de grande taille regardons la comparaison UMFPACK/GC-SSOR :

  • En 1d le GC est très compétitif cela est du au très faible nombre d’itération pour converger - avantage UMFPACK - Il est possible d’implémenter un pivot de Gauss spécifique pour une matrice tridiagonale (Algorithme TDMA) qui peut-être encore plus efficace.
  • Le 2d UMFPACK garde l’avantage, le GC met trop d’itérations le préconditionneur est ici peu efficace… Attention on ne compte pas le côut mémoire…
  • Pour le 3d cette fois la comparaison est en faveur du GC ! Et ceci à double titre : le GC fait mieux ou aussi bien à toutes les tailles de matrice et en plus il permet d’atteindre des dimensions non-envisageable avec UMFPACK.

Le Gradient Conjugué préconditionné SSOR est une méthode facile d’implémentation. Sont efficacité à $w_{opt}=\frac{2}{1+\sin(\pi h)}$ le rend compétitif en 3d. Quand serait-il avec un préconditionneur ILU, SPAI ou Multigrille…

Remarque : JULIA est un jeune langage qui quelquefois surprend par ces temps d’éxécution. Notez la réécriture de certaines fonctions afin d’obtenir de meilleurs performances. Enfin cette dernière courbe donnant l’avantage au gradient conjugué en 3d ne peut pas être obtenue sous MATLAB/OCTAVE/SCILAB…