Iterative Methods A X b Direct Methods Gauss

Iterative Methods A. X = b Ø Direct Methods : Gauss, LU, L LT , etc. Op(N 3) Ø Iterative methods : generates a sequences of vectors X(0), X(1), X(2), …X(n) such that: 1

Iterative Methods The term “Iterative Method” refers to a wide range of techniques that use successive approximations to obtain more accurate solutions to a linear system at each step. 2

Iterative Methods Østationary iterative method (Jacobi, Gauss-Seidel, SOR, etc): Iterative method that performs in each iteration the same operations on the current iteration vectors. . ØNonstationary iterative Method (Conjugate Gradient (CG)): Iterative method that has iteration-dependent coefficients. 3

Jacobi Method (Jacobi 1830) • The Jacobi method is based on solving for every variable locally with respect to the other variables; on iteration of the method corresponds to solving for every variable once. The resulting method is easy understand implement, but convergence is slow. 4

Jacobi Method (Jacobi 1830) To solve the nxn linear system A. x = b A = M + N M represent the diagonal part of A. N represent the strictly lower-triangular, and the strictly upper-triangular parts of A (M(-1) exists) 5

• A. x = b → M. x + N. x = b • Then x = -M-1. N. x + M-1. b Consider the iteration: x(0) : Choose the initial guess x(k+1) = -M-1. N. x(k) + M-1. b x(k+1) = Ω. x(k) + h With Ω = -M-1. N and h = M-1. b 6

• Theorem: For any x(0) the sequence {x(k)} (k = 0, 1, 3…) defined by: x(k+1) = Ω. x(k) + h For each K > 0 converges to the unique solution of x = Ω. x + h, and only if ρ(Ω ) < 1 • Definition: Diagonally dominant • Theorem: If A is stricly diagonally dominant, then for any choice of x(0) both the jacobi and Gauss. Seidel methods give sequences {x(k)} (k = 0, 1, 3…) that converge to the solution A. X=b 7

• In matrix terms: x(k+1) = Ω. x(k) + h • The ith iteration of the method defined: • We obtain 8

• Example: 9

10

Example 11

• The form x(k+1) = Ω. x(k) + h 12

• For an initial approximation let x(0)=(0, 0, 0, 0)T and generate x(1) by (Jacobi method): 13

14

• Results: 15

• For an initial approximation let x(0)=(0, 0, 0, 0)T and generate x(1) by (Gauss-Siedel method): 16

17

18

Gauss-Siedel method (Von Siedel 1830) • The Gauss Seidel method is like the Jacobi method, except that it uses updated values as soon as they available. In general, if the Jacobi method converges, the Gauss-Seidel method will converge faster than the Jacobi method, Though still relatively slowly. 19

Gauss-Siedel Method (Von Siedel 1830) To solve the nxn linear system A. x = b A = M + N 20

• A. x = b → M. x + N. x = b Then x = -M-1. N. x + M-1. b x(0) initial approximation x(k+1) = -M-1. N. x(k) + M-1. b x(k+1) = Ω. x(k) + h with Ω = -M-1. N and h = M-1. b 21

• In matrix terms: x(k+1) = Ω. x(k) + h • The ith iteration of the method defined: • We obtain 22

Stopping criterion • Since an iterative method computes successive approximations to the solution of a linear system, a practical test is needed to determine when to stop the iteration. Ideally this test would measure the distance of the last iterate to the true solution, but this is not possible. Instead, various other metrics are used, typically involving the residual. 23

• Stopping criterion: • Virtual Convergence: 24

(SOR: Successive Over Relaxation): • Successive Over Relaxation can be derived from the Gauss-Seidel method by introducing an extrapolation parameter w. For the optimal choice of w, SOR may converge faster than Gauss-Seidel by an order of magnitude. 25

• successive over relaxation (SOR): • If A is stricly diagonally dominant, and if 0< w <1 then for any choice of x(0) the SOR method give sequences {x(k)} (k = 0, 1, 3…) that converge to the solution A. X=b. • If A SDP(symmetric positive definite ), and 0 < w < 2. Then SOR method converge. 26

Choosing the value of ω • If ω = 1, the SOR method simplifies to the GS. • In general, it is not possible to compute in advance the value of ω that is optimal with respect to the rate of convergence of SOR. • Some heuristic estimate is used, such as ω = 2 – O(h) where h is the mesh spacing of the discretization of the underlying physical problem. • Theoretically optimal value (PDEs): 27

28

• The equations for the SOR method with ω =1. 25 29

30

31 Méthodes itératives stationnaires

Accélérateur Aiken • Aitken: une série géométrique basée sur les deux itérations précédentes. • La suite Xi est le vecteur solution donné par la méthode de Jacobi et de Gauss-Seidel. La suite Yi est le vecteur solution par Aitken. 32

Méthodes itératives stationnaires Ou nonstationnaires ? ? Les méthodes de projection qu’on verra par la suite sont plus générales et plus robustes que les méthodes stationnaires. 33

Méthodes itératives non-stationnaires Notations Kk(A, r 0) = eng{r 0, Ar 0 , …. . , Ak-1 r 0} espace de Krylov associé à A et r 0 x* est la solution exacte & le résidu r(i) = A x(i) – b Gradient conjugué : (Conjugate Gradient (CG)) La méthode du CG est une méthode efficace pour les matrices SDP 34

Théorie du (CG) x(i) élément de x(0) + Kk(A, r 0) E(x) = (x(i)-x*)T A (x(i)-x*) soit minimisé La résolution du système A x = b est remplacée par une minimisation : J(x) = ½ xt A x – xt b On cherche les itérés x(k+1) sous la forme : x(k+1) = x(k) + m(k) u(k) ) Où m(k) minimise E dans la direction u(k) choisie. 35

• Conjugate Gradient (CG): The conjugate gradient method derives its name from the fact that it generates a sequence of conjugate (or orthogonal) vectors. These vectors are the residuals of the iterates. They are also the gradients of a quadratic functional, the minimization of which is equivalent to solving the linear system. CG is an extremely effective method when the coefficient matrix is symmetric positive definite, since storage for only a limited number of vectors is required. 36

The CG Method Choose an initial guess x(0) to the solution x. r(0) = A x(0) – b u(0) = - r(0) a(0) = (r(0) , r(0)) Iterates : m(k) = a(k) _ (u(k) , A u(k)) x(k+1) = x(k) + m(k) u(k) r(k+1) = A x(k+1) – b = r(k) + m(k) A u(k) a(k+1) = (r(k+1) , r(k+1)) (End if a(k+1) = 0) l(k+1) = a(k+1) a(k) u(k+1) = - r(k+1) + l(k+1) u(k) k = 0, … , n-1 37

• Example: The linear system A. x = b 38

39

40

• Convergence Accurate predictions of the convergence of iterative methods are difficult to make. Spectral condition number: 41

Convergence du (CG) Ø Il est difficile de prédire exactement la convergence des méthodes itératives. Ø L’erreur se borne en fonct⁰ du conditionn. spectral k 2 de A Ø Si A est SDP, avec un préconditionneur M SDP, on peut démontrer que : ||x(i)-x*||A 2 ai || x(0) – x*||A Avec : Ø Si les valeurs propres extrémales de la matrice M-1 A sont très éloignées, On aura un phénomène qui s’appelle convergence superlinéaire. 42

Implémentation du (CG) La méthode du gradient conjugué implique à chaque itération : Ø Un produit matrice-vecteur de A par la direction de descente u Ø Trois opérations de mise à jour pour les vecteurs : Résidu, direction de descente et vecteur solution Ø Deux produits scalaires (u(k) , A u(k)) et (r(k+1) , r(k+1)) Lorsqu’on préconditionne le système par exemple par le préconditionneur SSOR, chaque produit matrice-vecteur sera transformé en une suite de produits matrice-vecteur : A * u sera D 1/2 (D + L)-1 A ( D + LT)-1 D 1/2 u : 43

Choix du vecteur de départ x 0 de CG : On peut affirmer que la méthode du gradient conjugué converge le mieux lorsque le vecteur de départ est fixé à zéro pour tous les systèmes linéaires Remarquons qu’on ne peut pas cependant généraliser ce résultat à tout type de matrices 44

Récapitulation des méthodes non-stationnaires : 45

Suite des méthodes non-stationnaires : 46

ST OC KA GE La méthode de stockage dépend de la méthode de résolution Ø La factorisation d’une matrice creuse détruit son caractère creux, sauf si elle est bande. Ø Une méthode itérative, n’utilise la matrice qu’à travers sa multiplication par des vecteurs alors on peut ne pas stocker que les coefficients non nuls de la matrice. Voire ne rien stocker. Soit A creuse de taille n et soit Nz le nombre total d’éléments non nuls 47

STOCKAGE (2) 48

Préconditionnement La convergence des méthodes itératives dépend du conditionnement spectral de A Avec : (k 2(A) = lmax(A) / lmin(A) ) A x = b M-1 A x = M-1 b k 2(A) > k 2(M-1 A) Compromis de coût entre construction & application du préconditionneur M à chaque itération et le gain dans la vitesse de convergence 49

Préconditionnement à gauche et à droite M-1 A n'est pas symétrique, même si M et A le sont. M = M 1 M 2 M 1 -1 AM 2 -1(M 2 x) = M 1 -1 b M 1 = le préconditionneur gauche et M 2 = le préconditionneur droit Et à la fin des itérations : x M 2 -1 xn Avec xn la solution finale calculée 50

Préconditionnement (1) M = D = M 1 M 2 = D 1/2 Préconditionneur de Jacobi D = Diag(A) M 1 -1 AM 2 -1(M 2 x) = M 1 -1 b (D-1/2 A D-1/2 ) (D 1/2 x) = (D-1/2 b) M = (D + L) D-1 ( D + L )T Ou Le préconditionneur SSOR M( ) = 1 2 - ( 1 D + L) ( 1 D)-1 ( 1 D + L)T A = D + LT Pour A Sym. M = (D + L) D-1/2 ( D + L )T 51

Préconditionneur SSOR Le truc d’Eisenstat’s <=> concevoir un nouveau système A’ u = b’ (D + L)-1 A ( D + LT)-1 ( D + LT) x = (D + L)-1 b Préconditionner A x = b avec M = (D + L) D-1 ( D + L )T Préconditionner A’ u = b’ avec le préconditionneur M = D-1 D 1/2 A’ D 1/2 ( D-1/2 u ) = D 1/2 b’ 52

Comparaison entre méthodes directes & méthodes itératives La méthode directe est commode pour les systèmes denses d’ordre supérieur, ainsi que pour les matrices bandes même d’ordre élevé. La méthode itérative est mieux adaptée aux autres matrices d’ordre très élevé et comportant de nombreux éléments nuls. Mais, lorsque le code cessera de tourner séquentiellement, et l’idée de la parallélisation s’impose, c’est la méthode itérative qui remporte la victoire. 53