Calcul de valeurs propres d'une matrice symétrique

 

Alors que la résolution d'un système linéaire passe par la détermination de l'unique solution à un système de la forme A.x = b, où A et b sont connus, on peut également être amené à résoudre des systèmes du type A.x = mx, où m est une constante qui constitue une inconnue supplémentaire du problème. En-dehors de la solution triviale x = 0, certaines valeurs de m (appelées valeurs propres) permettent de déterminer une famille de vecteurs x solution du système, appelés vecteurs propres.

La méthode la plus immédiate pour résoudre ce genre de systèmes consiste à résoudre le polynôme en m qui résulte de l'écriture du déterminant de (A-m.I). C'est pourtant loin d'être la solution la plus efficace dès que la taille de la matrice A augmente, d'autant plus que la détermination des zéros d'un polynôme de degré élevé présente également des difficultés. Là encore méthodes directes et itératives s'opposent, et on commencera cette fois par exposer ces dernières.

Méthode de la puissance

Permet de déterminer la plus grande valeur propre (en valeur absolue) de la matrice A et le vecteur propre associé.
Partant d'un vecteur initial arbitraire x⁽⁰⁾ dont une composante est unitaire, on construit un vecteur intermédiaire y⁽¹⁾ tel que A.x⁽⁰⁾ = y⁽¹⁾, puis on construit x⁽¹⁾ en "renormalisant" y(1) de manière à ce que la composante choisie soit à nouveau égale à l'unité. L'opération doit être répétée jusqu'à la convergence, en changeant éventuellement de composante de référence si celle-ci converge vers 0.
Le facteur de renormalisation converge vers la plus grande valeur propre et x converge vers le vecteur propre associé, sous les conditions que la plus grande valeur propre soit unique et que tous les vecteurs propres soient indépendants. La vitesse de convergence est proportionnelle au rapport entre les deux plus grandes valeurs propres (en valeur absolue).

Considérons par exemple le système [[3 1][1 3]] et partons du vecteur [0 1]. Le premier vecteur intermédiaire est [[3 1][1 3]] .[0 1] = [1 3], que l'on renormalise en [1/3 1]. On a ensuite [[3 1][1 3]] . [1/3 1] = [2 10/3], que l'on renormalise en [3/5 1]. L'étape suivante donne [[3 1][1 3]] . [3/5 1] = [14/5 18/5], renormalisé en [7/9 1], etc ... La série ainsi construite converge effectivement vers le vecteur propre [1 1], et le facteur de renormalisation (successivement égal à 3, 3.33, 3.6 ...) converge vers la valeur propre 4.

Méthode de la puissance inverse

L'idée est de déterminer la plus grande valeur propre (en v.a) de la matrice inverse A^-1, dont l'inverse est la plus petite valeur propre (toujours en v.a) de la matrice A.
Procédons comme pour la méthode de la puissance mais en remplaçant évidemment A par A^-1 lors de la construction de la suite de vecteurs intermédiaires :
A^-1.x^(k) = y^(k+1)
En multipliant cette équation par A on obtient :
A.A^-1x^(k) = A.y^(k+1)
que l'on peut écrire A.y^(k+1) = x^(k), système linéaire que l'on va résoudre à l'aide d'une décomposition LU. Une fois que y^(k+1) est déterminé, on procède comme précédemment à une renormalisation pour calculer les itérés suivants x^(k+1), qui convergeront vers le vecteur propre recherché ...

Méthode de la puissance avec décalage

Les valeurs propres d'une matrice A peuvent être décalées d'une valeur s en soustrayant s.x au problème A.x = m.x. On obtient (A-s.I).x = (m-s).x. Le décalage de A par s.I provoque un décalage d'une quantité s de la valeur propre m, mais les vecteurs propres restent inchangés. L'utilité de cette technique est triple :
- supposons qu'une matrice dont les deux valeurs propres, égales à 3 et 7, doivent être déterminées. En choisissant un décalage égal à 7, les deux valeurs propres de la matrice décalée deviennent -4 et 0.  L'application de la méthode de la puissance permet de déterminer la plus grande valeur propre (en v.a), à savoir -4. On en déduit par décalage inverse que la plus petite valeur propre de la matrice initiale est -4+7 = 3, et on obtient ainsi le même résultat qu'avec la puissance inverse ...
-  en choisissant un décalage inférieur à la plus grande vp, on peut de la même manière déterminer des valeurs propres intermédiaires.
- après plusieurs itérations de la méthode de la puissance, on peut utiliser un décalage égal à l'approximation obtenue de la plus grande valeur propre. La matrice décalée aura donc une plus petite valeur propre proche de l'unité, et on convergera plus rapidement vers la solution exacte en appliquant ensuite la méthode de la puissance inverse. 

Méthode directe

Les méthodes itératives ne permettent pas de résoudre des problèmes aux valeurs propres non-linéaires, c'est-à-dire définis par un problème de la forme A.X = B(m).x, où B(m) est une fonction non-linéaire de m. On doit alors passer par la résolution du polynôme caractéristique, mais ceci ne donne que les valeurs propres. On doit réinjecter ces dernières dans le système ou appliquer la méthode de la puissance inverse pour déterminer le vecteur propre associé.

Méthode QR

Contrairement aux méthodes présentées jusque ici, la méthode QR permet de déterminer simultanément toutes les valeurs propres (mais pas les vecteurs, encore une fois). Elle utilise le processus d'orthonormalisation de Gram-Schmidt pour transformer la matrice A en une matrice triangulaire similaire, dont les valeurs propres se retrouvent sur la diagonale.
Prenons l'exemple de la matrice [[3 4][4 3]].
Les vecteurs colonne de A sont a₁ = {3 4} et a₂ = {4 3}.
On crée d'abord le vecteur q₁ = a₁ / ||a₁|| = {3/5 4/5}, puis on construit le vecteur a'₂ = a₂ - (q₁^t.a₂)q₁ perpendiculaire à q₁. Ici, a'₂ = {28/25 -21/25}. La normalisation de a'₂ donne q₂ = {28/35 -21/35}.
La matrice Q de la décomposition QR est constituée des vecteurs q₁ et q₂, donc Q =  [[3/5 28/35][4/5 -21/35]]
La matrice R, pour sa part, est définie par r_ii = ||a'_i|| pour les termes diagonaux et par r_ij = q_i^t.a_j pour les termes extra-diagonaux, soit r₁₁ = 5, r₂₂ = 35/25 et r₁₂ = 24/5. On peut vérifier que A = Q.R ...
On utilise ensuite le fait que les produits QR et RQ donnent des matrices similaires. En effet, si A = QR, alors Q^-1A = Q^-1QR = R et Q^-1AQ = RQ.
A=QR et RQ ont donc des valeurs propres identiques, ce dont on peut ensuite profiter en appliquant de nouveau le processus ci-dessus à la matrice RQ. La matrice Q converge ainsi vers I, et les produits QR successifs convergent vers une matrice triangulaire ...

Autres méthodes

La méthode QR est généralement considérée comme étant la plus robuste et la plus performante pour déterminer l'ensemble des valeurs propres d'une matrice quelconque. Des méthodes plus efficaces permettent cependant de traiter les matrices "particulières". La méthode de Jacobi permet par exemple de transformer de manière itérative une matrice symétrique en une matrice diagonale. Les méthodes de Given et de Householder réduisent pour leur part de manière directe une matrice symétrique en une matrice tridiagonale, dont on pourra ensuite résoudre l'équation caractéristique par des techniques itératives.
La méthode de Householder peut cependant également être appliquée à une matrice non-symétrique, que la réduction sous forme de matrice de Hessenberg permettra de traiter ensuite par la méthode QR.

 

 

Equations et systèmes d'équations non-linéaires

Certaines méthodes de détermination des racines d'une équation non-linéaire (ou d'un système) passent par un encadrement de ces dernières, alors que d'autres conduisent à travailler sur un domaine non-borné. On compte les méthodes de la bissection et de la regula falsi au nombre des premières, et les méthodes du point fixe, de Newton, de la sécante au nombre des secondes.
Quelques principes généraux peuvent être énoncés :
- le choix d'une approximation initiale influe considérablement sur la rapidité de la convergence, voire sur la convergence elle-même
- les méthodes basées sur un domaine non-borné sont moins robustes car elles ne garantissent pas une convergence aussi systématique que les méthodes d'encadrement, mais quand elles convergent leur convergence est généralement plus rapide

 

Bissection

Si l'on a deux valeurs a et b telles que f(a).f(b) < 0, on calcule f((a+b)/2) et on modifie une des bornes de l'intervalle de recherche en fonction du résultat obtenu. Exemple simplissime : f(x) = (x-2). On part de a = 0 et b = 5. Comme f(a).f(b) < 0, on calcule f((a+b)/2). Puisque f(a).f((a+b)/2) <0, la racine est comprise entre a et (a+b)/2. Le processus doit être répété jusqu'à l'obtention de la précision voulue.
La quantité |b-a|/|a+b| est une approximation de l'erreur relative. Si L est la longueur de l'intervalle de recherche initial, la longueur après n itérations est évidemment de L/2ⁿ. On peut donc déterminer dès le début le nombre d'itérations nécessaire pour obtenir une erreur absolue donnée.
La bissection échoue dans les cas suivants :
- un nombre pair de racines dans l'intervalle de recherche initial.
- une racine en laquelle la dérivée première s'annule (ex : f(x) = x² en 0)

 

Regula falsi

Le choix qui consiste à calculer une nouvelle valeur de la fonction au point milieu de l'intervalle de recherche dans la méthode de la bissection n'est pas forcément le plus pertinent, en particulier dans le cas où la valeur en l'une des bornes de l'intervalle est beaucoup plus grande que la valeur en l'autre borne :  |f(a)| >> |f(b)|
On peut par exemple faire l'hypothèse que f est approchée par une fonction g linéaire entre les points a et b :
g'(x) = [f(b)-f(a)] / (b-a)
La racine c de g est donnée par :
c = b - f(b)/g'(x)
On remplace ensuite comme dans la bissection l'une des deux bornes de l'intervalle de recherche par c selon le signe de f(a).f(c) et de f(b).f(c) ...

 

Méthode du point fixe

On transforme f(x) = 0 en g(x) = x, et on utilise l'itération x_n+1 = g(x_n). On peut étudier la convergence de cette méthode en faisant un développement de l'erreur. Soit e_n = x_n - r l'erreur à l'itération n. On a

e_n+1 = g(r+e_n) - g(r) = g'(r)e_n + O(e_n²)

L'erreur à l'étape n+1 est donc directement proportionnelle à l'étape à l'itération n, et ne pourra diminuer que si |g'(r)| < 1. La convergence est dans ce cas linéaire, sauf cas particuliers ou les dérivées n-ièmes g⁽ⁿ⁾(r) sont nulles, auquel cas la convergence est d'ordre n+1 (on rappelle que l'ordre de convergence d'une méthode est p si l'on a |e_n+1| = C |e_n|^p + O(e_n^p+1) ...).
Remarque : puisque e₁ ~ g'(r)e₀ et que e₂ ~ g'(r)e₁, on a e₂/e₁ ~ e₁/e₀. On peut donc gagner un ordre de convergence avec ce qui constitue l'extrapolation d'Aitken (algo de Steffenson) : r ~ (x₂x₀ - x₁²)/(x₂-2x₁+x₀), ou plutôt, par souci de stabilité : r ~ x₀ - (x₁-x₀)²/(x₂-2x₁+x₀).

 

Méthode de Newton

x_n+1 = x_n - f(x_n)/f '(x_n)

C'est une méthode de point fixe particulière, avec g(x) = x - f(x)/f '(x) ... Comme pour celle-ci, la méthode converge si l'estimation initiale est suffisamment proche de la racine, mais la convergence est ici quadratique (dans le cas de racines simples) car on a en outre g'(r) = 0. Dans le cas de racines multiples, on peut récupérer la convergence quadratique en remplaçant f par u = f / f '. Toute racine multiple de f devient racine simple de u.
(Rappel : la convergence est dite quadratique - resp. d'ordre n - quand l'erreur commise lors d'une itération est proportionnelle au carré - resp. à la puissance nième - de l'erreur commise à l'itération précédente) 

Lorsque f est une fonction "simple", f ' peut être déterminée analytiquement. En revanche, dans le cas général d'une fonction non-linéaire la dérivée doit être approximée, par exemple selon
f '(x) = (f(x+e) - f(x)) / e
ce qui double le nombre d'évaluations de la fonction. Des erreurs d'arrondi peuvent en outre apparaître si e est trop petit; alors que la convergence peut être ralentie si e est trop grand.
Lorsque l'on a à résoudre un système d'équations, on adapte la méthode précédente en remplaçant f '(x_n) par la matrice jacobienne Df(x) définie par Df(x)_ij = df_i(x)/dx_j ( il s'agit ici de "d rond"). Compte tenu du coût représenté par l'évaluation de la dérivée (ou du jacobien), coût encore plus sensible avec un système d'équations que dans le cas d'une équation unique, on peut choisir de conserver la même valeur de f '(x), ou de Df(x)_ij, pendant plusieurs itérations. Cette méthode de Newton approchée voit disparaître la convergence quadratique. On peut également (méthodes de quasi-Newton) ne réaliser un calcul exact du jacobien que toutes les p itérations, mais modifier (sans le recalculer exactement) ce dernier à chaque itération.

 

Méthode de la sécante

On peut simplifier la méthode de Newton en remplaçant le calcul de dérivée par l'approximation f '(x_n) = (f(x_n) - f(x_n-1)) / (x_n -x_n-1). Il faut fournir au départ deux valeurs initiales : c'est un algorithme à deux pas.
On peut montrer que la convergence quadratique est perdue, mais qu'elle demeure plus que linéaire. On a en fait |e_n+1| ~ C.|e_n|^1,618. On reconnaît dans l'exposant le nombre d'or (sqrt(5)+1)/2 ...
Il a été montré (Jeeves) que si l'effort calculatoire nécessaire pour évaluer f ' est inférieur à 0.43 fois l'effort nécessaire à l'évaluation de f alors la méthode de Newton est plus performante. Dans le cas contraire la méthode de la sécante sera préférée.

 

Méthode de la corde

C'est une méthode de Newton approchée dans laquelle on évite le calcul systématique de f '(x_n) en remplaçant cette dernière par f '(x₀) à chaque itération. La convergence est alors linéaire.
Dans le cas d'un système, cette méthode permet de ne réaliser qu'une fois le coûteux calcul de jacobien, et si l'on peut faire une décomposition LU ou LL^t de D, on n'aura que deux systèmes triangulaires à résoudre à chaque itération. Reste à savoir si le gain en temps compense la perte en vitesse de convergence ...

 

Méthode de Muller

Alors que la méthode de la sécante consistait à approcher f par une fonction linéaire au voisinage de la racine, la méthode de Muller consiste à l'approcher localement par une fonction quadratique. Trois valeurs initiales sont donc nécessaires ... L'ordre de convergence de cette méthode est de 1.84, ce qui la situe donc entre les méthodes de la sécante (1.62) et de Newton (2). On estime généralement que le gain apporté par rapport à la méthode de la sécante ne justifie pas le surcoût calculatoire.

 

Détermination des zéros d'un polynôme

Les méthodes présentées ci-dessus peuvent être utilisées pour déterminer les racines d'un polynôme, mais le traitement des racines multiples peut se révéler problématique. On peut par exemple adapter la méthode de Newton, soit en tenant compte de la multiplicité m de la racine pour modifier l'itération de la manière suivante :

x_n+1 = x_n - m.f(x_n)/f '(x_n)

soit en étudiant la fonction modifiée u(x) = f(x) / f '(x). On montre aisément que la racine unique de u est également la racine multiple de f. L'application de la méthode de Newton à u(x) passe cependant par le calcul de f ''(x).

Lorsque les racines sont susceptibles d'être complexes, les méthodes basées sur les techniques d'encadrement sont mises en défaut, simplement parce que le signe de f(x) ne change pas forcément lorsque la racine est complexe. Les méthodes de Newton, de Muller ou de la sécante devront être adaptées à l'utilisation d'une arithmétique complexe.