Optimisation quadratique

On s’intéresse dans ce chapitre à la résolution numérique des systèmes d’équations linéaires par des méthodes itératives finies ; nous verrons plus loin ce que l’on entend par ces deux qualificatifs. Pour certaines méthodes, on considérera leur extension à la résolution de systèmes d’équations non linéaires. Pour d’autres, cette extension se fera dans d’autres chapitres. On cherche donc un point x∗ ∈ R n solution du système linéaire en x suivant Ax = b, (8.1) où A est une matrice d’ordre n donnée et b est un vecteur donné dans R n. Nous supposerons toujours que A est inversible. Dans ce cas, la solution unique de (8.1) s’écrit x∗ = A −1 b, où A−1 est la matrice inverse de A. Pour obtenir x∗, on ne calcule jamais A−1 . Cela n’est ni nécessaire, ni opportun dès que n dépasse quelques unités (le temps de calcul requis est alors trop important), ni stable numériquement [313 ; 2002, section 14.2]. Rappelons que la résolution d’un système linéaire d’ordre n non structuré par factorisation de la matrice A demande de l’ordre de O(n 3 ) opérations. L’algorithme le plus simple et le plus utilisé pour ce faire est l’élimination gaussienne, équivalente à la factorisation gaussienne de la matrice A, avec pivotage partiel. Cet algorithme requiert 2n 3/3 + O(n 2 ) opérations flottantes. Les algorithmes présentés dans cette section ont aussi une complexité opérationnelle en O(n 3 ). Nous distinguerons deux cas : celui où A est symétrique définie positive et celui où A est non symétrique. Si la matrice A est symétrique, la solution du système (8.1) est aussi le point stationnaire de la fonction quadratique f(x) = 1 2 x TAx − b Tx. (8.2) Si, de plus, A est définie positive, x∗ réalise le minimum de f. On peut dans ce cas obtenir des algorithmes de résolution du système (8.1) à partir d’algorithmes de minimisation de f. L’algorithme du gradient conjugué est de ce type (section 8.2). Si A n’est pas symétrique, on peut remplacer la résolution du système (8.1) par la résolution de l’équation normale équivalente A TAx = A T b. (8.3) Comme la matrice de ce système est symétrique définie positive, on peut songer à utiliser les techniques mentionnées ci-dessus. Cependant, le système (8.3) peut être beaucoup moins bien conditionné que le système original (8.1), si bien qu’il est souvent préférable d’utiliser des méthodes s’attaquant directement au système (8.1), sans faire appel à des algorithmes de minimisation de fonction. L’idée est alors de générer des approximations de x∗ qui sont, dans un sens à préciser, les meilleures approximations sur des sous-espaces affines de dimension croissante. On peut en effet espérer que si l’on a une « bonne » approximation de x∗ sur le sous-espace affine x1+Kp, Kp étant un sous-espace vectoriel de dimension p, il ne sera pas trop difficile d’obtenir une « bonne » approximation sur un sous-espace affine x1 + Kp+1, de dimension p + 1, contenant x1 + Kp. Dans l’algorithme du résidu minimal (section 8.3), les sous-espaces Kp sont obtenus par un procédé particulier : ce sont des sous-espaces de Krylov (section 8.1). D’ailleurs, nous verrons que l’algorithme du gradient conjugué génère des itérés qui sont aussi dans un tel sous-espace de Krylov, si bien que ces deux algorithmes font partie de la même famille, celle dite des méthodes de Krylov. Les méthodes de Krylov ont la propriété de trouver la solution du système (8.1) en un nombre fini d’étapes (en artihmétique exacte). Elles s’apparentent sur ce point aux méthodes directes de résolution, fondées sur la factorisation de la matrice A : factorisation gaussienne, factorisation QR… En pratique cependant, la présence d’erreurs d’arrondi dans les calculs empêche les méthodes de Krylov de trouver la solution en un nombre fini d’étapes, dès que la dimension du système linéaire dépasse quelques dizaines, si bien qu’on les range souvent dans la famille des méthodes itératives. Ces dernières sont intéressantes pour plusieurs raisons. D’abord, comme elles procèdent par améliorations successives, on peut s’arrêter lorsque la précision est jugée satisfaisante. Ceci peut se produire bien avant que la solution ne soit trouvée. Une solution « acceptable » peut donc être obtenue à un faible coût. Ensuite, elles permettent de tirer parti d’une bonne approximation initiale, ce qui est souvent le cas lorsqu’on doit résoudre une succession d’équations linéaires provenant de la linéarisation d’une même équation non linéaire. Elles ont aussi des inconvénients, par exemple, de ne pouvoir exploiter la « creusité » (caractère creux) éventuelle de A que par des produits matrice-vecteur plus rapides et d’avoir, en pratique, besoin d’un bon préconditionneur pour converger en un nombre raisonnable d’itérations. Notons pour terminer qu’il existe aussi des méthodes itératives ne trouvant pas la solution en un nombre fini d’itérations : méthodes de Jacobi, de Gauss-Seidel, de relaxation… Les premières ont été introduites au XIXe siècle pour calculer plus rapidement (à la main, bien sûr) des solutions approchées, alors que la résolution par élimination directe était trop fastidieuse (voir la lettre de Gauss en épigraphe de ce chapitre). Elles sont cependant généralement considérées comme moins efficaces que les méthodes de Krylov [601]. Pour une introduction à ces méthodes, on pourra consulter les livres de Varga [605 ; 1962], de Ciarlet [126 ; 1982] et de Hackbusch [298 ; 1994]. Connaissances supposées. L’introduction de l’algorithme du gradient conjugué se fait en utilisant le procédé d’orthogonalisation de Gram-Schmidt (section B.1.1). 

Sous-espaces de Krylov 

On appelle sous-espace de Krylov d’ordre p, associé à une matrice A d’ordre n et à un vecteur r ∈ R n, le sous-espace vectoriel de R n engnedré par les vecteurs r, Ar, 8.1. Sous-espaces de Krylov 343 A2 r, . . . , Ap−1 r. On le note Kp ≡ Kp(A, r) ≡ vect{r, Ar, . . . , Ap−1 r}. La proposition suivante montre que la dimension de Kp vaut p, tant que p reste inférieur ou égal à un certain indice s, à partir duquel sa dimension ne bouge plus (et vaut s). On peut avoir s < n. On dit que le sous-espace de Krylov Ks est saturé et que s = s(A, r) est l’indice de saturation des sous-espaces de Krylov associés à A et à r. La proposition montre aussi que l’indice de saturation est atteint dès que A−1 r est « capturé » par les sous-espaces de Krylov. Proposition 8.1 Si A est inversible, alors il existe un indice s > 1 tel que K1 $ · · · $ Ks = Kp, ∀p > s. L’indice s est caractérisé par le fait que A −1 r ∈ Ks \ Ks−1. 

Proposition 8.2

Si A est inversible, alors pour tout r ∈ R n s(A, r) 6 β, où β est de degré du polynôme minimal unitaire annihilant A. Cette majoration est optimale : on peut trouver un vecteur r pour lequel on a s(A, r) = β.Les deux algorithmes itératifs de résolution du système linéaire (8.1) que nous décrivons dans ce chapitre, l’algorithme du gradient conjugué (section 8.2) et l’algorithme GMRES (section 8.3), sont des méthodes de Krylov, ce qui veut dire que l’itété xk appartient à un certain sous-espace de Krylov Kk(A, r). Dans les deux cas, le vecteur r est le résidu r1 = b − Ax1 évalué en l’itéré initial x1. La raison pour laquelle ce choix de r convient est expliqué au début de la section 8.3.1. L’algorithme du gradient conjugué est adapté aux systèmes linéaires dans lesquels la matrice A est symétrique définie positive (ou semi-définie positive), alors que l’algorithme GMRES est plus général (et plus difficile à mettre en œuvre) puisqu’il s’affranchit de toute hypothèse sur A. 

Notion de directions conjuguées

 Pour une matrice A symétrique définie positive donnée, la notion de conjugaison équivaut à la notion d’orthogonalité pour le produit scalaire associé à A. De façon plus précise, on a les définitions suivantes. On dit que u et v ∈ R n sont conjuguées si u TAv = 0. Plus généralement, on dit que u1, . . . up ∈ R n sont conjuguées si u T i Auj = 0, ∀i 6= j. Enfin, on dit que u est conjuguée par rapport à un sous-espace vectoriel E ⊆ R n si u TAv = 0, ∀v ∈ E. Si v1, . . . , vp sont des vecteurs linéairement indépendants dans R n, on pourra leur associer p directions conjuguées u1, . . . , up en utilisant le procédé du Gram-Schmidt (section B.1.1) avec le produit scalaire hu, vi = u TAv, dans lequel on se passe des étapes de normalisation superflues (étapes 1 et 2.3). Il suffit de prendre, pour k = 1, . . . , p

Algorithme des directions conjuguées

Une méthode de directions conjuguées pour la minimisation d’une fonction quadratique f est une méthode qui à chaque itération prend xk ∈ x1+Ek de façon à minimiser f sur le sous-espace affine x1+Ek. Le lien avec la conjugaison des directions est donné par la proposition suivante. On note gk = ∇f(xk) = Axk − b le gradient de f en xk pour le produit scalaire euclidien et yk = gk+1 − gk.