Valeurs propres

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• expliquer pourquoi le polynôme caractéristique n'est pas un algorithme de calcul des valeurs propres, en invoquant à la fois l'impossibilité d'Abel–Ruffini et le conditionnement catastrophique des racines d'un polynôme donné par ses coefficients;
• localiser le spectre d'une matrice dans le plan complexe avec les disques de Gershgorin, démontrer le théorème et utiliser ses deux raffinements (composante connexe de $k$ disques, version transposée);
• mettre en œuvre la méthode de la puissance, démontrer sa convergence vers le vecteur propre dominant et prédire son coût à partir du rapport $|{\lambda }_2/{\lambda }_1|$;
• améliorer une estimation de valeur propre par le quotient de Rayleigh et justifier sa précision quadratique sur une matrice symétrique;
• atteindre n'importe quelle valeur propre par la puissance inverse avec décalage, en reconnaissant qu'elle coûte une résolution de système par itération — donc une factorisation $LU$ faite une fois pour toutes;
• énoncer le principe de l'algorithme $QR$ et le théorème de Bauer–Fike, et dire pourquoi une matrice non symétrique peut avoir des valeurs propres d'une sensibilité extravagante.

Pourquoi le polynôme caractéristique n'est pas un algorithme

Le problème de ce chapitre est celui-ci: étant donné $A\in {\mathbb{R}}^{n\times n}$, trouver les couples $(\lambda ,\mathbf{v})$ avec $\mathbf{v}\ne 0$ tels que

L'algèbre linéaire répond en une ligne: $\lambda$ est valeur propre si et seulement si $A-\lambda I$ est singulière, c'est-à-dire si et seulement si

Le polynôme caractéristique $p_A$ est de degré $n$; ses racines sont les valeurs propres; le problème est «résolu». Un cours d'analyse numérique doit commencer par expliquer pourquoi cette réponse, parfaitement correcte, ne donne aucun algorithme utilisable. Il y a deux raisons, et elles sont indépendantes l'une de l'autre.

Première raison: il n'existe pas de formule

Depuis Abel (1824) et Ruffini, on sait qu'il n'existe pas de formule par radicaux donnant les racines du polynôme général de degré $n\ge 5$. Galois a ensuite expliqué pourquoi: le groupe symétrique $S_n$ n'est pas résoluble pour $n\ge 5$. Ce n'est pas une lacune de notre ingéniosité, c'est un théorème d'impossibilité.

La conséquence pour nous est structurelle. Toute méthode de calcul des valeurs propres d'une matrice $n\times n$ avec $n\ge 5$ est nécessairement itérative: elle produit une suite qui converge vers le spectre, et on l'arrête quand on est satisfait. Il n'existe pas, et il ne peut pas exister, d'analogue de l'élimination de Gauss — laquelle, elle, résout $A\mathbf{x}=\mathbf{b}$ en un nombre fini d'opérations connu d'avance. C'est une différence de nature entre le chapitre 3 et celui-ci, et elle vaut la peine d'être méditée: la résolution d'un système linéaire est un problème fini, le calcul d'un spectre ne l'est pas.

Seconde raison: les racines d'un polynôme sont mal conditionnées

On pourrait croire qu'il suffit alors de former $p_A$, puis d'appliquer une méthode itérative de recherche de racines — Newton, par exemple, dont le chapitre 2 a montré qu'elle double les chiffres corrects à chaque pas. C'est là que se trouve le vrai obstacle, et il est numérique, pas algébrique: passer par les coefficients d'un polynôme détruit l'information.

Rappelons la mesure du conditionnement d'un problème introduite au chapitre 1: le facteur par lequel une erreur relative sur la donnée est multipliée dans le résultat. Soit $p$ un polynôme de degré $n$ écrit $p(x)={\sum }_{k=0}^n{a}_k{x}^k$ et soit une racine de . Perturbons le seul coefficient en ; la racine devient . En dérivant implicitement par rapport à , on obtient , d'où

Le facteur ${\kappa }_k(r)$ est le conditionnement de la racine $r$ vis-à-vis du coefficient $a_k$. Rien n'empêche $a_k{r}^k$ d'être gigantesque pendant que est modeste — et c'est exactement ce qui arrive dès que le polynôme a plusieurs racines de tailles voisines. Wilkinson a rendu ce phénomène célèbre en 1963 sur le polynôme de degré 20 dont les racines sont ; on l'observe déjà très bien au degré 5.

Le cas d'une racine multiple est encore plus brutal, et il se dit en une ligne. Pour $p(x)=(x-1{)}^5$, une perturbation $\epsilon$ du terme constant donne $(x-1{)}^5=\epsilon$, donc $x-1={\epsilon }^{1/5}$: les racines se déplacent comme la de la perturbation. Avec , les cinq racines quittent d'exactement — un déplacement fois plus grand que la cause, et elles deviennent complexes, réparties sur un petit cercle. Le conditionnement d'une racine multiple est infini: (5.3) le voit d'ailleurs directement, puisque annule le dénominateur.

Le problème aux valeurs propres vu comme un point fixe

Puisqu'il faut itérer, autant reformuler (5.1) sous la forme qui appelle une itération. Un vecteur propre n'est pas un vecteur, c'est une direction: si $\mathbf{v}$ convient, $\alpha \mathbf{v}$ convient aussi pour tout $\alpha \ne 0$. La bonne inconnue est donc un point de la sphère unité, et la bonne application est l'application normalisée.

Un vecteur propre est donc un point fixe de $\Phi$ — au signe près, qui est le prix de la normalisation par une norme. Le chapitre 2 nous a appris quoi faire d'un point fixe: itérer, et regarder si l'itération contracte. C'est littéralement la méthode de la puissance, et le taux de contraction sera $|{\lambda }_2/{\lambda }_1|$.

Deux conventions pour toute la suite. On note ${\lambda }_1,\dots ,{\lambda }_n$ les valeurs propres de $A$, rangées par module décroissant,

et $\rho (A)=|{\lambda }_1|$ le rayon spectral, déjà utilisé au chapitre 4 pour caractériser la convergence des méthodes itératives. Les normes vectorielles et matricielles, le rayon spectral et les propriétés des matrices symétriques sont rappelés dans l'annexe A; nous les utilisons ici sans les redémontrer.

Localiser le spectre: les disques de Gershgorin

Avant de calculer une valeur propre, il est souvent bien plus utile de savoir où elle se trouve. Le théorème de Gershgorin (1931) donne cette information gratuitement — sans une seule itération, à la simple lecture des coefficients. Il sert à choisir un décalage, à prouver qu'une matrice est inversible, à borner un rayon spectral, à certifier qu'un schéma itératif converge.

Démonstration. Soit $\lambda$ une valeur propre de $A$ et $\mathbf{v}\ne 0$ un vecteur propre associé. Puisque $\mathbf{v}$ est non nul, on peut choisir un indice $k$ où son module est maximal:

Écrivons la $k$-ième ligne de l'égalité $A\mathbf{v}=\lambda \mathbf{v}$:

On isole le terme diagonal en le faisant passer à droite:

Prenons les modules et appliquons l'inégalité triangulaire au membre de gauche:

où la seconde majoration utilise $|v_j|\le |v_k|$ pour tout $j$ — c'est précisément pour pouvoir écrire cela qu'on a choisi $k$ à la composante maximale. Comme $|v_k|>0$, on peut diviser par :

c'est-à-dire $\lambda \in D_k$. La valeur propre appartient donc bien à la réunion des disques. $\square$

La démonstration dit plus que l'énoncé, et il faut le remarquer: le disque qui contient $\lambda$ est celui dont l'indice correspond à la plus grande composante du vecteur propre. C'est ce qui rend le résultat pratique — une valeur propre dont le vecteur propre est concentré sur la coordonnée $k$ est piégée dans $D_k$, et ce disque est petit dès que la ligne $k$ est presque diagonale.

Les deux raffinements

Démonstration (esquisse). On relie $A$ à sa diagonale par le chemin continu

de sorte que $A(0)=D$ et $A(1)=A$. Les disques de $A(t)$ ont les mêmes centres que ceux de $A$ et des rayons $t{r}_i$, donc ils sont contenus dans ceux de $A$ pour tout : la réunion des disques considérés reste, pour tout , disjointe de celle des autres. Or les valeurs propres d'une matrice dépendent de ses coefficients (elles sont les racines d'un polynôme dont les coefficients sont polynomiaux en ceux de la matrice, et les racines d'un polynôme dépendent continûment de ses coefficients). En , la réunion contient exactement les coefficients diagonaux correspondants, donc valeurs propres. Quand croît de à , aucune valeur propre ne peut franchir la frontière de la réunion sans quitter, à cet instant, tous les disques de — ce que le théorème 5.1 interdit. Le compte est donc conservé, et vaut en .

Ce raffinement admet une démonstration complète par le théorème des résidus ou par la continuité des valeurs propres au sens de la multiplicité algébrique; nous nous en tenons à l'argument de déformation, qui contient toute l'idée. Le cas $k=1$ est celui qu'on utilise le plus: un disque isolé contient exactement une valeur propre, et comme le spectre d'une matrice réelle est stable par conjugaison, cette valeur propre est alors nécessairement réelle — sans quoi sa conjuguée, qui est aussi valeur propre, serait dans le même disque, qui en contiendrait deux.

Démonstration. Une matrice et sa transposée ont le même polynôme caractéristique, puisque $\det (A^{\mathsf{T}}-\lambda I)=\det ((A-\lambda I{)}^{\mathsf{T}})=\det (A-\lambda I)$: elles ont donc exactement le même spectre. Or les disques de lignes de sont, par définition, les disques de colonnes de . Le théorème 5.1 appliqué à donne le résultat, et l'intersection (5.7) parce que les deux inclusions valent simultanément.

C'est un raffinement gratuit, et l'exemple 5.2 a montré qu'il peut diviser un rayon par trois. On retiendra la règle: calculer les deux jeux de disques coûte $2n^2$ additions et ne se refuse pas.

Une conséquence immédiate mérite d'être isolée, car elle relie ce chapitre aux deux précédents. Si $A$ est à diagonale strictement dominante, c'est-à-dire si $|a_{ii}|>r_i$ pour tout $i$, alors $0$ n'appartient à aucun disque $D_i$ — car . Donc n'est pas valeur propre, et . C'est le théorème de Hadamard, obtenu ici en une ligne; c'est aussi l'hypothèse sous laquelle le chapitre 4 démontre la convergence de la méthode de Jacobi, et ce n'est pas une coïncidence: la matrice d'itération de Jacobi est , dont les disques de Gershgorin sont centrés en et de rayon , ce qui donne directement.

La méthode de la puissance

L'idée et l'algorithme

Reprenons le point fixe (5.4). Si $A$ est diagonalisable, ses vecteurs propres ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ forment une base, et tout vecteur de départ s'y décompose:

Appliquons $A$ un grand nombre de fois. Chaque composante est multipliée par sa propre valeur propre:

Si $|{\lambda }_1|$ domine strictement les autres modules, chaque rapport ${\lambda }_i/{\lambda }_1$ est de module $<1$ et sa puissance $k$-ième tend vers zéro: la direction de $A^k{\mathbf{x}}^{(0)}$ s'aligne sur . Tout le contenu de la méthode est dans le crochet de (5.8), et le terme qui s'éteint le moins vite est celui d'indice 2.

Reste le problème d'échelle: le facteur ${\lambda }_1^k$ fait déborder ou disparaître les nombres. Avec ${\lambda }_1=3,2143$ et $k=800$, ${\lambda }_1^k$ dépasse et vaut en double précision; avec , il vaut bien avant. , ce qui ne change rien à la direction et retire le facteur parasite. Nous normalisons par la composante de plus grand module, convention qui présente deux avantages: le vecteur garde une composante égale à exactement — donc aucune ambiguïté de signe — et le facteur par lequel on a divisé converge vers .

def puissance(A, x, n):
    """n iterations de la methode de la puissance; renvoie (mu, x)."""
    for _ in range(n):
        y = [sum(A[i][j] * x[j] for j in range(len(x))) for i in range(len(A))]
        j = max(range(len(y)), key=lambda i: abs(y[i]))
        mu 

Le coût par itération est celui d'un produit matrice-vecteur: $n^2$ multiplications et $n^2-n$ additions, soit $2n^2$ opérations, plus $n$ divisions pour la normalisation. C'est le point fort de la méthode: $O(n^2)$ par itération au lieu de , et seulement si la matrice est creuse avec un nombre borné de coefficients par ligne. C'est la raison pour laquelle la méthode de la puissance, vieille de cent ans, reste la brique des algorithmes de classement de très grands graphes, où se compte en milliards et où aucune factorisation n'est envisageable.

La convergence et son taux

Démonstration. Écrivons la décomposition de l'itéré non normalisé, comme en (5.8):

La mise en facteur est licite parce que $c_1\ne 0$ (hypothèse 2) et ${\lambda }_1\ne 0$ (hypothèse 1, qui force $|{\lambda }_1|>0$ dès que ). Majorons le reste par l'inégalité triangulaire:

la dernière majoration utilisant $|{\lambda }_i|\le |{\lambda }_2|$ pour tout $i\ge 2$, c'est-à-dire exactement le rangement par module décroissant. L'hypothèse 1 donne $q=|{\lambda }_2/{\lambda }_1|<1$, donc géométriquement de raison .

Passons à la direction. Le vecteur normalisé s'écrit

où ${\sigma }_k$ est un scalaire de module $1$ — c'est lui qui porte l'ambiguïté de signe: si ${\lambda }_1<0$, ${\sigma }_k$ alterne, et ${\mathbf{x}}^{(k)}$ bascule d'un côté à l'autre de la même droite. La distance à la droite engendrée par ne voit pas , et il suffit d'évaluer

Posons $\mathbf{u}={\mathbf{v}}_1+{\mathbf{w}}_k$. Pour deux vecteurs non nuls, l'identité

et l'inégalité $|\Vert \mathbf{u}\Vert -\Vert {\mathbf{v}}_1\Vert |\le \Vert \mathbf{u}-{\mathbf{v}}_1\Vert =\Vert {\mathbf{w}}_k\Vert$ donnent

Pour $k$ assez grand, $\Vert {\mathbf{w}}_k\Vert \le \frac{1}{2}\Vert {\mathbf{v}}_1\Vert$, donc $\Vert \mathbf{u}\Vert \ge \frac{1}{2}\Vert {\mathbf{v}}_1\Vert$ et le majorant vaut au plus . Avec (5.11) on obtient (5.9), la constante absorbant le facteur numérique. La distance décroît donc au moins comme .

Pour la valeur propre, écrivons la $j$-ième composante de $A{\mathbf{x}}^{(k)}$ à partir de (5.10). En posant ${\mathbf{u}}_k={\mathbf{v}}_1+{\mathbf{w}}_k$, on a et donc

d'où, en prenant le rapport des $j$-ièmes composantes,

Le numérateur est majoré par $(\Vert A\Vert +|{\lambda }_1|)\Vert {\mathbf{w}}_k\Vert$, donc en $O(q^k)$; le dénominateur tend vers $(v_1{)}_j\ne 0$, car est l'indice de la composante maximale, qui est non nulle à la limite. Ainsi .

Le message de la figure 5.3 est économique. Le nombre d'itérations est $\ln (\text{tol})/\ln r$: il est modeste tant que $r$ reste loin de 1 et il explose ensuite. À $r=0,2$ il faut une dizaine d'itérations; à $r=0,9$, cent vingt et une; à $r=0,99$, plus de treize cents. Autrement dit, la méthode de la puissance est excellente quand une valeur propre se détache nettement, et inutilisable quand deux valeurs propres sont voisines — ce qui est, malheureusement, le cas le plus fréquent en pratique. Toute la suite du chapitre consiste à fabriquer artificiellement un grand écart là où il n'y en avait pas.

Le quotient de Rayleigh

Le rapport ${\mu }_k$ de deux composantes utilise une seule coordonnée de l'itéré et jette les $n-1$ autres. On peut faire beaucoup mieux, et le résultat est l'une des plus jolies estimations de l'analyse numérique.

L'interprétation géométrique est immédiate, et c'est celle que l'explorateur ci-dessous met en scène: pour $\mathbf{x}$ unitaire, $\rho (\mathbf{x})={\mathbf{x}}^{\mathsf{T}}(A\mathbf{x})$ est la longueur signée de la projection de $A\mathbf{x}$ sur $\mathbf{x}$. C'est «de combien $A$ étire $\mathbf{x}$ dans sa propre direction», en ignorant ce qui est parti de côté. Une autre lecture, tout aussi utile: $\rho (\mathbf{x})$ est la valeur de $\alpha$ qui minimise , c'est-à-dire pour le vecteur donné. C'est cette optimalité qui explique le théorème suivant.

Démonstration. Le théorème spectral (algèbre linéaire; rappelé à l'annexe A) fournit une base orthonormée de vecteurs propres, dans laquelle on écrit $\mathbf{x}={\sum }_i{\alpha }_i{\mathbf{q}}_i$ avec ${\sum }_i{\alpha }_i^2=1$ puisque est unitaire et la base orthonormée. Alors

Le quotient de Rayleigh est donc une moyenne pondérée des valeurs propres, de poids ${\alpha }_i^2\ge 0$ de somme 1: il est nécessairement compris entre ${\min }_i{\lambda }_i$ et ${\max }_i{\lambda }_i$, ce qui est la dernière assertion. Retranchons :

le terme $i=1$ s'annulant identiquement. D'où

Il reste à relier ${\sum }_{i\ge 2}{\alpha }_i^2$ à la distance au vecteur propre. Quitte à changer $\mathbf{x}$ en $-\mathbf{x}$, on peut supposer ${\alpha }_1\ge 0$, et alors

l'inégalité $2-2{\alpha }_1\ge 1-{\alpha }_1^2=(1-{\alpha }_1)(1+{\alpha }_1)$ étant vraie car . En reportant, on obtient (5.13).

Le point décisif est l'annulation du terme d'indice 1 dans la somme: la dérivée du quotient de Rayleigh est nulle en un vecteur propre. Un vecteur propre est un point critique de $\rho$, et l'erreur au voisinage d'un point critique est quadratique — c'est le même mécanisme qui fait qu'un extremum d'une fonction lisse est bien mieux localisé que sa position ne l'est.

La puissance inverse, et le décalage qui change tout

La méthode de la puissance ne sait atteindre que ${\lambda }_1$, et seulement si elle se détache. Deux transformations élémentaires suffisent à lever les deux limitations.

Inverser

Si $A$ est inversible, les valeurs propres de $A^{-1}$ sont les $1/{\lambda }_i$, avec les mêmes vecteurs propres — il suffit de multiplier $A\mathbf{v}=\lambda \mathbf{v}$ par ${\lambda }^{-1}{A}^{-1}$. La plus grande en module est : appliquer la méthode de la puissance à fait converger vers la valeur propre de de , au taux .

Décaler

Les valeurs propres de $A-\sigma I$ sont les ${\lambda }_i-\sigma$, toujours avec les mêmes vecteurs propres. En combinant, les valeurs propres de $(A-\sigma I{)}^{-1}$ sont les

et la plus grande en module correspond à la valeur propre ${\lambda }_i$ la plus proche de $\sigma$. Voilà la clé: en choisissant $\sigma$, on désigne la valeur propre que l'on veut, et on la rend dominante autant qu'on le souhaite.

Le taux de convergence se lit sur (5.14). Si ${\lambda }_p$ est la valeur propre la plus proche de $\sigma$ et ${\lambda }_s$ la suivante,

Ce rapport tend vers $0$ quand $\sigma \to {\lambda }_p$: plus le décalage est bon, plus la convergence est rapide, et sans borne. C'est le renversement complet de la situation de la figure 5.3.

Le coût, et pourquoi le chapitre 3 vient avant

Chaque itération de (5.15) demande la résolution d'un système linéaire au lieu d'un produit matrice-vecteur: $O(n^3)$ au lieu de $O(n^2)$, si l'on procède naïvement. La parade est celle du chapitre 3: la matrice $A-\sigma I$ ne change pas d'une itération à l'autre tant que $\sigma$ est fixé. On la factorise donc une seule fois, en $\frac{2}{3}{n}^3$ opérations avec pivot partiel, puis chaque itération se réduit à une descente et une remontée triangulaires, opérations. Le surcoût par rapport à la méthode de la puissance directe est nul en régime, et le coût fixe est celui d'une factorisation.

C'est exactement pour cela que le chapitre 3 précède celui-ci: sans $PA=LU$, la puissance inverse serait inabordable. Et c'est aussi la raison pour laquelle on ne calcule jamais $(A-\sigma I{)}^{-1}$ explicitement — trois fois plus cher, et numériquement moins stable, pour une information dont on n'utilisera que le produit par des vecteurs.

Le choix du décalage est donc l'enjeu. Trois sources, dans l'ordre de raffinement: les disques de Gershgorin, qui coûtent $2n^2$ additions; une valeur approchée connue d'avance (une fréquence propre attendue, un mode physique); et surtout, la mise à jour dynamique $\sigma \leftarrow \rho ({\mathbf{x}}^{(k)})$, dite itération du quotient de Rayleigh. Cette variante recalcule le décalage à chaque pas avec la meilleure estimation disponible; elle converge cubiquement pour une matrice symétrique — trois fois plus de chiffres corrects à chaque itération — au prix d'une nouvelle factorisation à chaque pas. C'est un compromis à peser: $\frac{2}{3}{n}^3$ par itération, mais trois ou quatre itérations en tout.

La déflation

Une fois ${\lambda }_1$ et ${\mathbf{v}}_1$ obtenus, comment atteindre ${\lambda }_2$ sans recommencer? L'idée de la déflation est de modifier $A$ de façon à remplacer ${\lambda }_1$ par en laissant les autres valeurs propres intactes; la méthode de la puissance, relancée, converge alors vers .

La vérification est immédiate sur ${\mathbf{v}}_1$: le terme retranché appliqué à ${\mathbf{v}}_1$ donne ${\lambda }_1{\mathbf{v}}_1({\mathbf{w}}_1^{\mathsf{T}}{\mathbf{v}}_1)/({\mathbf{w}}_1^{\mathsf{T}}{\mathbf{v}}_1)={\lambda }_1{\mathbf{v}}_1$, donc . Pour les autres, on utilise le fait que les vecteurs propres à gauche et à droite associés à des valeurs propres sont orthogonaux: si avec , alors et .

Sur la matrice témoin, avec ${\lambda }_1=3,2143$, ${\mathbf{v}}_1=(-0,272782,0,668755,1{)}^{\mathsf{T}}$ et le vecteur propre à gauche , la formule (5.18) donne

dont le spectre, recalculé, vaut $\{-0,675130870567;0,460811127189;-1,3\cdot 10^{-16}\}$: les deux valeurs propres restantes sont retrouvées à la précision machine et la troisième est numériquement nulle, comme annoncé.

Le principe de l'algorithme QR

Toutes les méthodes précédentes calculent une valeur propre à la fois. La méthode de référence pour obtenir le spectre complet d'une matrice dense est l'algorithme $QR$, publié indépendamment par Francis et par Kublanovskaya au début des années 1960, et toujours au cœur des bibliothèques d'algèbre linéaire soixante ans après. Nous en donnons le principe et les faits d'usage; l'analyse de sa convergence est admise, avec sa raison.

L'algorithme repose sur une observation d'algèbre linéaire: deux matrices semblables, $A$ et $P^{-1}AP$, ont le même polynôme caractéristique, donc le même spectre. On peut donc transformer $A$ autant qu'on veut par similitude sans rien perdre — et l'idée est de la transformer jusqu'à ce que le spectre se lise sur la diagonale.

Le point essentiel tient en une ligne: puisque $R_k=Q_k^{\mathsf{T}}{A}_k$, la nouvelle matrice s'écrit

c'est-à-dire que $A_{k+1}$ est semblable à $A_k$ par une transformation orthogonale. Toutes les matrices de la suite ont donc exactement le même spectre que $A$, et la transformation est orthogonale, donc parfaitement conditionnée au sens de la norme euclidienne. Sous des hypothèses de séparation des modules, la suite $A_k$ converge vers une forme triangulaire supérieure (triangulaire par blocs si est réelle avec des valeurs propres complexes), dont la diagonale porte les valeurs propres.

Trois faits d'usage, que toute bibliothèque met en œuvre.

1. Réduction préalable à la forme de Hessenberg. Telle quelle, chaque itération coûte $\frac{4}{3}{n}^3$: inabordable. On commence donc par rendre $A$ semblable, par des réflexions de Householder, à une matrice de Hessenberg (nulle sous la première sous-diagonale) en $\frac{10}{3}{n}^3$ opérations. Cette forme est par l'itération , et une itération sur une Hessenberg ne coûte plus que . Si est symétrique, la forme de Hessenberg est tridiagonale et une itération coûte .

Ce que le cas symétrique apporte

Quand $A$ est symétrique réelle, le théorème spectral change complètement la situation: il existe une matrice orthogonale $Q$ et une matrice diagonale réelle $\Lambda$ telles que $A=Q\Lambda Q^{\mathsf{T}}$. Les conséquences sont toutes bonnes, et il vaut la peine de les lister:

• toutes les valeurs propres sont réelles — plus de paires complexes, plus de blocs $2\times 2$, plus de problème de dominance pour la méthode de la puissance;
• les vecteurs propres sont orthogonaux entre eux, donc la base propre est parfaitement conditionnée (${\kappa }_2(Q)=1$);
• le quotient de Rayleigh est quadratiquement précis (théorème 5.5) et encadré par les valeurs propres extrêmes;
• la forme de Hessenberg est tridiagonale, donc l'itération $QR$ coûte $O(n)$ par pas;
• enfin, et c'est le plus important, les valeurs propres sont parfaitement conditionnées, au sens de la section suivante.

C'est pourquoi les bibliothèques numériques distinguent systématiquement le cas symétrique, avec des routines séparées: ce n'est pas une optimisation, c'est un autre problème.

Le conditionnement d'une valeur propre

Reste la question qui ouvrait le chapitre, posée cette fois sur la matrice et non sur le polynôme: si l'on perturbe $A$ en $A+\delta A$, de combien le spectre bouge-t-il?

Démonstration. Si $\mu$ est l'une des ${\lambda }_i$, l'inégalité est triviale. Supposons donc $\mu \notin \mathrm{sp}(A)$, de sorte que $\Lambda -\mu I$ est inversible. Comme $\mu$ est valeur propre de $A+\delta A$, la matrice est singulière; en la conjuguant par ,

est singulière elle aussi. En factorisant à gauche le facteur inversible $\Lambda -\mu I$,

est singulière, donc le crochet l'est. Or une matrice de la forme $I+M$ singulière impose $\Vert M\Vert \ge 1$ pour toute norme induite — sinon la série de Neumann ${\sum }_k(-M{)}^k$ convergerait et fournirait un inverse. Donc

La matrice $(\Lambda -\mu I{)}^{-1}$ est diagonale, de coefficients $1/({\lambda }_i-\mu )$; sa norme 2 est le plus grand module de ceux-ci, soit $1/{\min }_i|{\lambda }_i-\mu |$. En reportant et en réarrangeant,

ce qui est (5.21). $\square$

Deux lectures, opposées.

Si $A$ est symétrique, le théorème spectral donne $V=Q$ orthogonale, donc ${\kappa }_2(Q)=1$, et (5.21) se lit

Le déplacement du spectre ne dépasse pas la perturbation elle-même: les valeurs propres d'une matrice symétrique sont parfaitement conditionnées, sans aucune hypothèse de séparation. Une perturbation de $10^{-16}$ déplace le spectre d'au plus $10^{-16}$. (Le théorème de Weyl dit même mieux: la $i$-ième valeur propre bouge d'au plus $\Vert \delta A{\Vert }_2$, valeur propre par valeur propre.)

Si $A$ n'est pas symétrique, ${\kappa }_2(V)$ peut être arbitrairement grand — et il l'est précisément quand deux vecteurs propres sont presque parallèles, c'est-à-dire quand la matrice est presque défective.

Synthèse

• Le polynôme caractéristique n'est pas un algorithme. Abel–Ruffini interdit une formule pour $n\ge 5$, donc toute méthode est itérative; et surtout, les racines d'un polynôme donné par ses coefficients sont souvent atrocement mal conditionnées — facteur $160$ au degré 5, $9\cdot 10^4$ au degré 10, infini pour une racine multiple. La traduction utile va dans l'autre sens: on calcule les racines d'un polynôme par l'algorithme $QR$ appliqué à sa matrice compagnon.
• Les disques de Gershgorin localisent le spectre pour $2n^2$ additions. Chaque valeur propre est dans la réunion des disques centrés sur la diagonale et de rayon la somme des modules hors diagonale. Un disque isolé contient exactement une valeur propre, donc réelle si est réelle; la version par colonnes donne une seconde localisation, à intersecter avec la première. Le théorème ne dit pas que chaque disque contient une valeur propre.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 5 (approximation des valeurs propres et des vecteurs propres).
• Trefethen, L. N. et Bau, D., Numerical Linear Algebra, SIAM, Philadelphie, leçons 24 à 29 (méthode de la puissance, itération du quotient de Rayleigh, algorithme $QR$ avec décalages).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chapitre sur le calcul des valeurs propres.
• Golub, G. H. et Van Loan, C. F., Matrix Computations, Johns Hopkins University Press, Baltimore, chap. 7 et 8 (le problème aux valeurs propres non symétrique et symétrique, et l'analyse complète de l'algorithme $QR$ admise ici).
• Wilkinson, J. H., The Algebraic Eigenvalue Problem, Clarendon Press, Oxford (l'ouvrage où le conditionnement des racines d'un polynôme est établi et illustré sur le polynôme de degré 20 qui porte son nom).
• Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 9.