On découvre ici différentes techniques de calcul des éléments propres d'une matrice donnée, ainsi que des applications associées. On suppose que le lecteur est familier avec les notions de diagonalisation.
Ce TP s'inspire du livre Maths et Maple de Jean-Michel Ferrard : il condense le chapitre 6, le transpose à python, et propose des travaux pratiques en plus.
Pour obtenir les valeurs propres, importez le paquetage linalg, et utilisez la fonction eigvals :
>>> from numpy import * >>> A = matrix((1,2,3,4,5,6,7,8,1)).reshape((3,3)) >>> eigvals(A) array([ 12.45416474, -0.37976219, -5.07440255])
On peut obtenir les vecteurs propres ainsi :
>>> eig(A)
(array([ 12.45416474, -0.37976219, -5.07440255]),
matrix([ [-0.29373774, -0.73967882, -0.29720654],
[-0.69005397, 0.66500848, -0.39870229],
[-0.66147083, -0.10314536, 0.86758559] ]))
Le retour est constitué du tableau des valeurs propres, suivi de la matrice des vecteurs propres.
Donnons un premier exemple de diagonalisation matricielle, dans le cas d'une matrice diagonalisable. On commence par importer ce qu'il faut, et définir la matrice à diagonaliser...
>>> from numpy import * >>> from numpy.linalg import * >>> A = matrix((1,2,3,4,5,6,7,8,9)).reshape((3,3)) >>> set_printoptions(suppress = True)
On rappelle que la fonction eig renvoie le tableau des valeurs propres, et la matrice des vecteurs propres.
>>> eig(A)
(array([ 16.11684397, -1.11684397, -0. ]),
matrix([[-0.23197069, -0.78583024, 0.40824829],
[-0.52532209, -0.08675134, -0.81649658],
[-0.8186735 , 0.61232756, 0.40824829]]))
$A$ possède 3 valeurs propres distinctes, et est de taille $3 \times 3$, donc est diagonalisable.
L'écriture diagonale de $A$ correspond à l'insertion des valeurs propres dans la diagonale d'une matrice nulle. Cela se fait grâce à la fonction diag:
>>> D=diag(eig(A)[0])
>>> D
array([[ 16.11684397, 0. , 0. ],
[ 0. , -1.11684397, 0. ],
[ 0. , 0. , -0. ]])
Quant à la matrice de passage P de la base canonique vers la base où $A$ est diagonale, elle correspond à la matrice des vecteurs propres :
>>> P=eig(A)[1]
On vérifie que l'on a bien $A = P D P^{-1}$ :
>>> P*D*inv(P)
matrix([[ 1., 2., 3.],
[ 4., 5., 6.],
[ 7., 8., 9.]])
>>> A
matrix([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
>>> D
array([[ 16.11684397, 0. , 0. ],
[ 0. , -1.11684397, 0. ],
[ 0. , 0. , -0. ]])
La diagonalisation est réalisée !
Reste à obtenir le polynôme caractéristique :
>>> from numpy import * >>> M=matrix((1,2,3,4)).reshape((2,2)) >>> poly(M) array([ 1., -5., -2.])
C'est-à-dire que le polynôme caractéristique est $\chi_M(\lambda) = \lambda^2-5\lambda-2$.
Pour obtenir les valeurs propres, on utilise la méthode eigenvals :
>>> from sympy import * >>> A = Matrix((-5,2,4,-4,
10,-3,-6,8,
11,-4,-6,8,
22,-8,-14,17)).reshape(4,4) >>> A [-5, 2, 4, -4] [10, -3, -6, 8] [11, -4, -6, 8] [22, -8, -14, 17] >>> vp=A.eigenvals() >>> vp {1: 2, 2: 1, -1: 1}
Une fois encore, le calcul est ici exact. Les valeurs derrière les : correspondent aux multiplicités des valeurs propres.
Le retour de la fonction eigenvals est un dictionnaire. Pour récupérer ses clés, qui correspondent aux valeurs propres à proprement parlé :
>>> vp.keys() [1, 2, -1]
Pour les vecteurs propres associés, la méthode est eigenvects :
>>> A.eigenvects()
Cependant, dans la version de sympy que je possède, cette méthode a un problème. Je vous propose un petit correctif :
>>> def eigenvects2(self, **flags):
... if flags.has_key('multiple'):
... del flags['multiple']
... out, vlist = [], self.eigenvals(**flags)
... for r, k in vlist.iteritems():
... tmp = self - eye(self.lines)*r
... basis = tmp.nullspace()
... if len(basis) != k:
... basis = tmp.nullspace(simplified=True)
... out.append((r, k, basis))
... return out
...
Pour l'utiliser :
>>> VP = eigenvects2(A) >>> len(VP) 3
Il y a donc trois valeurs propres différentes...
>>> VP[0] (1, 2, [[1] [1] [1] [0], [0] [2] [0] [1]])
Donc 1 est de multiplicité 2, et possède pour vecteurs propres linéairement indépendants $(1,1,1,0)^T$ et $(0,2,0,1)^T$. De même,
>>> VP[1] (2, 1, [[ 0] [ 1] [1/2] [ 1]])
Donc 2 est de multiplicité 1, et possède pour vecteur propre $\left(0,1,\frac{1}{2},1\right)^T$. etc.
Reste à montrer comment obtenir le polynôme caractéristique :
>>> P = A.charpoly('x')
>>> P
Poly(x**3 - 15*x**2 - 18*x, x)
On reprend l'exemple précédent de diagonalisation, en l'enrichissant : on va tester si la matrice considérée est diagonalisable, avant de réaliser sa diagonalisation proprement dite.
On rappelle quelles sont les conditions pour qu'une matrice $A$ puisse s'écrire sous la forme $P D P^{-1}$ où $D$ est diagonale, c'est-à-dire pour qu'il existe une base dans laquelle $A$ a une forme diagonale.
Pour cela, il faut que :
On rappelle que les sous-espaces propres sont les sous-espaces vectoriels engendrés par les vecteurs propres associés à chacune des valeurs propres (il y a un sous-espace propre pour chacune des valeurs propres distinctes).
On montre, dans ce qui suit, deux applications classiques de la diagonalisation : le calcul de la puissance nième d'une matrice, et la résolution d'un système différentiel...
Pour calculer $M^n$, on évite de faire $M \times M \times ... \times M$, $n$ fois, vu le coût d'un calcul matriciel.
On préfère la technique suivante, quand elle est possible : on diagonalise $M$, on élève la réduite diagonale $D$ à la puissance $n$, puis on utilise la relation
$M^n = P D^n P^{-1}$
où $P$ est la matrice de passage de la diagonalisation.
Voici un exemple : considérons la matrice $M$ suivante, dont on souhaite obtenir la puissance $M^{7}$...
$\left(\begin{array}{ccc}1 & 4 & 3 \\ 1 & -2 & 1 \\ 0 & 1 & 3 \end{array}\right)$
On calcule sa matrice diagonale, et sa matrice de passage :
>>> from numpy import *
>>> from numpy.linalg import *
>>> M = matrix((1,4,3,1,-2,1,0,1,3)).reshape((3,3))
>>> M
matrix([[ 1, 4, 3],
[ 1, -2, 1],
[ 0, 1, 3]])
>>> D,P=diag(eig(M)[0]),eig(M)[1]
>>> D
array([[-3.03404184, 0. , 0. ],
[ 0. , 1.48653522, 0. ],
[ 0. , 0. , 3.54750661]])
>>> P
matrix([[ 0.6505674 , 0.96068514, 0.87265468],
[-0.74922933, 0.23164273, 0.23451873],
[ 0.12416708, -0.15305459, 0.42833955]])
Alors la matrice $M^{7}$ sera égale à $P D^{7} P^{-1}$, où $D^{7}$ s'obtient en élevant chaque nombre sur la diagonale à la puissance 7 :
>>> D7 = zeros((3,3))
>>> D7[0,0]=D[0,0]**7
>>> D7[1,1]=D[1,1]**7
>>> D7[2,2]=D[2,2]**7
>>> D7
array([[-2366.74219863, 0. , 0. ],
[ 0. , 16.04081493, 0. ],
[ 0. , 0. , 7070.7013837 ]])
D'où la puissance septième de $M$ :
>>> P*D7*inv(P)
matrix([[ 656., 4392., 10664.],
[ 768., -1208., 2968.],
[ 440., 1648., 5272.]])
$\left(\begin{array}{ccc} 0 & 1 & 1 \\-3 & 4 & 3 \\ -1 & 1 & 2 \end{array}\right)$
$\left(\begin{array}{cc} 1 & 4 \\-1 & -2 \end{array}\right)$
On peut utiliser la diagonalisation de $M$ pour résoudre le système différentiel $X' = MX$.
En effet, si $M$ est diagonalisable, alors $ X' = MX$ est équivalent à $X' = P D P^{-1} X$. Donc $P^{-1} X' = D P^{-1} X$, et le système est alors équivalent à
$Y' = D Y$.
avec $Y = P^{-1} X$. Or, le nouveau système est d'une forme bien plus sympathique :
$\frac{d}{dt} \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \left( \begin{array}{cccc} d_1 & 0 & ... & 0 \\ 0 & d_2 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & d_n \end{array} \right) \times \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right)$
C'est-à-dire qu'il correspond au système
$\left\{ \begin{array}{ccc}y_1^' & = & d_1 y_1 \\ y_2^' & = & d_2 y_2 \\ \vdots & &\vdots \\ y_n^' & = & d_n y_n \end{array}\right.$
Ce système admet pour solution
$\left\{ \begin{array}{ccc}y_1 & = & e^{d_1t} u_1 \\ y_2 & = & e^{d_2t} u_2 \\ \vdots & &\vdots \\ y_n & = & e^{d_nt} u_n \end{array}\right.$
où $U = (u_1, ..., u_n)^T$ est le vecteur des constantes d'intégration, déterminables à l'aide des conditions initiales. Cette solution donne, sous forme matricielle :
$\left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \left( \begin{array}{cccc} e^{d_1t} & 0 & ... & 0 \\ 0 & e^{d_2t} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & e^{d_nt }\end{array} \right) \times \left( \begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_n \end{array} \right) = \left( \begin{array}{cccc} \sum_{k=0}^{+\infty} \frac{d_1^kt}{k!} & 0 & ... & 0 \\ 0 & \sum_{k=0}^{+\infty} \frac{d_2^kt}{k!} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & \sum_{k=0}^{+\infty} \frac{d_n^kt}{k!} \end{array} \right) \times \left( \begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_n \end{array} \right)$
On reconnaît alors une exponentielle de matrice...
$Y = \sum_{k=0}^{+\infty} \frac{1}{k!} \left( \begin{array}{cccc} d_1t & 0 & ... & 0 \\ 0 & d_2t & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & d_nt \end{array} \right)^k \times \left( \begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_n \end{array} \right) = e^{Dt} U$
Finalement, en se rappelant que $Y = P^{-1} X$...
Théorème : La solution du système différentiel $X' = MX$, où la matrice $M$ est diagonalisable ($M = PDP^{-1}$), est $X = P e^{Dt} U$, où $U = P^{-1}X(0)$.
On souhaite résoudre le système différentiel :
$\left\{ \begin{array}{lll} x_1^' & = & x_1 + 4 x_2 \\ x_2^' & = & x_1 - 2 x_2 \end{array} \right.$
avec la condition initiale $x_1(0) = 1, x_2(0) = 1$.
Définissons ce système en python :
>>> from sympy import * >>> M = Matrix((1,4,1,-2)).reshape(2,2) >>> X0 = Matrix((1,1)).reshape(2,1) >>> X0 [1] [1]
Puis, recherchons les éléments propres de $M$
>>> M.eigenvects() [(2, 1, [[4] [1]]), (-3, 1, [[-1] [ 1]])] >>> D = Matrix((2,0,0,-3)).reshape(2,2) >>> D [2, 0] [0, -3] >>> P = Matrix((4,-1,1,1)).reshape(2,2) >>> P [4, -1] [1, 1]
On peut maintenant calculer $e^{Dt}$ :
>>> t = Symbol('t')
>>> d = Matrix((exp(2*t),0,0,exp(-3*t))).reshape(2,2)
>>> d
[exp(2*t), 0]
[ 0, exp(-3*t)]
Puis $U$ :
>>> U = P.inv()*X0 >>> U [2/5] [3/5]
On peut donc obtenir la solution $X(t)$ :
>>> X = P*d*U >>> X [-3*exp(-3*t)/5 + 8*exp(2*t)/5] [ 3*exp(-3*t)/5 + 2*exp(2*t)/5]
On peut alors, par exemple, regarder quelle sera la solution du système aux temps $t=1$ et $t=7,5$ :
>>> X.subs(t,1) [-3*exp(-3)/5 + 8*exp(2)/5] [ 3*exp(-3)/5 + 2*exp(2)/5] >>> X.subs(t,1).evalf() [11.7926175172683] [2.98549468059298] >>> X.subs(t,7.5) [8*exp(15.0000000000000)/5 - 3*exp(-22.5000000000000)/5] [2*exp(15.0000000000000)/5 + 3*exp(-22.5000000000000)/5] >>> X.subs(t,7.5).evalf() [5230427.79595538] [1307606.94898884]
$\left\{ \begin{array}{lll} x_1^' & = & x_2 + x_3 \\ x_2^' & = & -3 x_1 + 4 x_2 + 3 x_3 \\ x_3^' & = & -x_1+x_2+2 x_3\end{array} \right.$
Les décompositions LU, QR et de Choleski peuvent servir à approcher les valeurs propres, voire les vecteurs propres...
Soit $A_1$ une matrice inversible. Alors on sait qu'il existe un couple $(Q_1,R_1)$ de matrices telles que $A_1 = Q_1 \times R_1$, où :
Soit alors $A_2 = R_1 \times Q_1$. On peut montrer que cette matrice est semblable à $A_1$ : elles ont les mêmes valeurs propres. On applique alors à $A_2$ une décomposition QR : $A_2 = Q_2 \times R_2$.
Soit alors $A_3 = R_2 \times Q_2$, etc.
Imaginons continuer ainsi le procédé. Alors on peut aboutir, sous certaines hypothèses, à trouver les valeurs propres de $A$.
Plus précisément, on peut montrer que :
Théorème :
Si :
* les valeurs propres de A sont de module deux à deux distincts,
* une matrice de passage diagonalisant A possède une décomposition
LU,
Alors la suite de matrices :
* $A_0 = A$,
* $A_{i+1} = R_i \times Q_i$, où $(Q_i,R_i)$ est le couple de la
décomposition QR de $A_i$
converge vers une matrice triangulaire, dont la diagonale est
constituée des valeurs propres de A.
Bien évidemment, les hypothèses du théorème précédent sont impossibles à vérifier, puisqu'on cherche précisément quelles sont les valeurs propres de $A$. On se contente donc d'agir de façon empirique.
Cette méthode peut facilement se mettre en œuvre, par exemple pour trouver les valeurs propres de la matrice
$\left( \begin{array}{ccc} -4 & 4 & -2 \\ 1 & -4 & 5 \\ -4 & 8 & 2 \end{array}\right)$
>>> from sympy import * >>> M = Matrix((-4,4,-2,1,-4,5,-4,8,2)).reshape(3,3) >>> for k in range(20): ... (Q,R) = M.QRdecomposition() ... M = R.evalf()*Q.evalf() ... >>> M [ -10.0000974898348, 2.13808816953320, 4.01150932531371] [ -5.21174585283222e-5, -1.99998604282076, 1.15476205049695] [-0.000361064881559241, 9.66940332900396e-5, 6.00008353265553]
On tend donc à se rapprocher d'une matrice triangulaire, dont la diagonale est approximativement égale à -10, -2, 6. Vérifions que l'on a bien cela pour valeurs propres...
>>> from numpy import * >>> from numpy.linalg import * >>> M = matrix((-4,4,-2,1,-4,5,-4,8,2)).reshape((3,3)) >>> eigvals(M) array([-10., -2., 6.])
C'est effectivement le cas.
On constate les points suivants :
$\left(\begin{array}{cccc} 337 & 124 & -101 & -11 \\ 528 & 453 & 2 & -178 \\ 1227 & 718 & -199 & -209 \\ 2037 & 1258 & -299 & -279 \end{array}\right)$
$\left(\begin{array}{ccc} 1 & 2 & 3 \\ 3 & 2 & 1 \\ 2 & 3 & 1 \end{array}\right)$
Dans le cas général,
on devra donc trouver mieux.
Cependant, cette méthode s'applique bien si la matrice est de Hessenberg (i.e. presque triangulaire supérieure, sauf les coefficients immédiatement sous la diagonale), ou tridiagonale symétrique.
Cela peut sembler franchement restrictif. Cependant, on peut montrer que toute matrice réelle (respectivement symétrique réelle) peut être transformée en matrice de Hessenberg (resp. tridiagonale symétrique), par des méthodes de type Householder, ou Givens.
Dans le même esprit, comme la rapidité de la convergence est fonction de la bonne séparation des modules des valeurs propres, on peut améliorer la méthode QR en effectuant un décalage de la matrice. On ne s'intéresse alors plus aux $A_i$, mais aux $A_i - \beta I_n$.
Nous ne détaillerons pas ici ces points. Voir le livre de Ferrard pour plus d'informations.
On peut, comme ci-dessus, utiliser la décomposition LU pour déterminer les valeurs propres d'une matrice $A$.
On rappelle que si, pour tout entier k compris entre 1 et la taille de $A$, la sous-matrice formée des $k$ premières lignes et des $k$ premières colonnes de $A$ est inversible, alors $A$ peut s'écrire sous la forme d'un produit matriciel $L \times U$, où
La méthode de déterminisation des valeurs propres est alors la suivante :
Sous certaines conditions non triviales, la suite des matrices $L_k$ tend vers la matrice identité, quand la suite $U_k$ tend vers une matrice triangulaire supérieure limite, dont les coefficients diagonaux sont les valeurs propres de $A$.
Cette technique peut s'adapter au cas où $A$ n'est pas LU-décomposable, en utilisant une certaine matrice de permutation.
$\left(\begin{array}{ccc} 4 & 4 & -2 \\ 1 & -4 & 5 \\ -4 & 8 & 2 \end{array}\right)$
On peut enfin utiliser la décomposition de Choleski, dernière décomposition matricielle classique, pour déterminer les valeurs propres d'une matrice.
Nous n'entrerons pas plus dans les détails.
Dans le TP sur les décompositions matricielles, on a vu une méthode de Givens permettant de décomposer une matrice carrée réelle inversible $M$ sous la forme $QR$.
Chacune de ces matrices de rotation est responsable de l'annulation d'un coefficient sous-diagonal de $M$, la progression s'effectuant :
On arrive donc à une écriture $GM = R$, où $G$ est le produit des matrices de rotation, et $R$ triangulaire supérieure à diagoale positive. Puis, à $M = QR$, où $Q = G^T = G^{-1}$.
On pourrait penser adapter cette méthode au calcul des valeurs propres et vecteurs propres, en espérant que :
Hélas, le deuxième point n'est pas vérifié, et la méthode ci-dessus ne convient pas.
Cependant, on peut quand même utiliser les rotations de Givens dans la recherche des valeurs propres de $M$, si $M$ est symétrique, ce qui sera supposé dans toute cette section.
Jacobi repose sur le fait que si $G$ est une matrice de rotation et $M$ symétrique, alors $M_1 = GMG^T$ est une matrice semblable à $M$ : elles ont les mêmes valeurs propres. Comme $M_1$ reste symétrique, on peut continuer le procédé avec $M_1$.
Le but est alors d'utiliser des matrices de rotation $G_k$, pour que la matrice $M_p$ issue du produit $G_p \times G_{p-1} \times ... \times G_1 \times M \times G_1^T \times ... \times G_{p-1}^T \times G_p^T$ est telle que ses valeurs propres soient bien plus faciles à déterminer en pratique que celles de $M$ (et, on le sait, $M$ et $M_p$ ont les mêmes valeurs propres).
On prendra pour $G_k$ des matrices de Givens bien choisies, de sorte que les 0 introduits conduisent à un produit $G_p \times G_{p-1} \times ... \times G_1 \times M \times G_1^T \times ... \times G_{p-1}^T \times G_p^T$ égal à une matrice diagonale : les valeurs propres de $M$ seront les coefficients de cette diagonale !
Dans le détail, à chaque nouvelle multiplication par une matrice de Givens,
Cela n'est pas grave : la méthode de Jacobi est ainsi faite que la somme des carrés des coefficients non diagonaux décroît strictement à chaque itération : on se rapproche donc, à chaque fois, un peu plus d'une matrice diagonale...
On rappelle la forme générale des matrices de rotation de Givens...
>>> from sympy import *
>>> def G(n,i,j,x):
... M = eye(n)
... M[i,i] = M[j,j] = cos(x)
... M[i,j] = -sin(x)
... M[j,i] = sin(x)
... return M
...
>>> x=Symbol('x')
>>> G(6,2,4,x)
[1, 0, 0, 0, 0, 0]
[0, 1, 0, 0, 0, 0]
[0, 0, cos(x), 0, -sin(x), 0]
[0, 0, 0, 1, 0, 0]
[0, 0, sin(x), 0, cos(x), 0]
[0, 0, 0, 0, 0, 1]
A chaque étape, on attaque le coefficient non diagonal $m_0$ de plus grande valeur absolue dans $M$, qui est en position $(i_0,j_0)$, en multipliant par une matrice de Givens $G$ choisie de manière à annuler $m_0$.
Cette matrice de Givens bien choisie $G$ est $G(n,i_0,j_0,x)$, où
On remplace alors la matrice $M$ par $G M G^T$, et on recommence.
Le procédé est reproduit :
Les valeurs propres (approchées) seront les éléments diagonaux de la dernière matrice obtenue.
Chaque coefficient annulé est en général désannulé par les étapes suivantes, mais l'important est qu'à chaque étape la matrice $M_{k+1}$ est plus proche d'une matrice diagonale que ne l'était $M_k$ : la méthode de Jacobi converge.
Réalisez la méthode de Jacobi.
Pour trouver les valeurs propres d'une matrice carrée de manière approchée, il existe donc plusieurs méthodes, plus ou moins efficaces, plus ou moins faciles à utiliser.
Certaines ne s'appliquent qu'à des matrices d'un type un peu particulier (symétriques réelles, etc.) Les méthodes de puissances itérées, couramment employées, ont l'avantage de s'appliquer à des matrices quelconques. Elles permettent de trouver à la fois des valeurs propres et des vecteurs propres.
Leur principal défaut est de supposer que les valeurs propres sont dans une configuration favorable, ce que l'on ne peut savoir à l'avance.
La méthode des puissances itérées directes permet d'approcher la valeur propre de plus grand module (encore faut-il qu'il y en ait une).
On peut alors envisager une méthode de déflation, pour trouver toutes les valeurs propres (encore faut-il que leurs modules soient deux à deux différents).
Notons, cependant, que cette méthode est source d'erreurs d'arrondis.
La méthode des puissances itérées inverses permet de trouver la valeur propre de plus petit module (sauf, éventuellement, 0).
Une variante de cette méthode permet alors d'approcher une valeur propre quelconque, à condition d'en connaître une estimation initiale.
Considérons une matrice diagonalisable, dont les valeurs propres ont des modules distincts deux à deux.
On part d'un vecteur quelconque, que l'on multiplie un certain nombre de fois par la matrice $M$ dont on cherche les éléments propres.
Alors, au bout d'un certain nombre d'itérations, le vecteur résultant tend vers un vecteur propre associé à la plus grande valeur propre (en module)...
Illustrons cette méthode sur la matrice
$M = \left(\begin{array}{ccc} -7 & -3 & 0 \\ 54 & 8 & -36 \\ 3 & -3 & -10 \end{array} \right)$
>>> from numpy import * >>> M = matrix((-7,-3,0,54,8,-36, 3, -3,-10)).reshape((3,3)) >>> from numpy.linalg import * >>> eigvals(M) array([ 2., -1., -10.])
Les valeurs propres à trouver sont donc 2, -1 et -10. Tirons un vecteur aléatoirement
>>> from random import randint
>>> U = matrix((0,0,0)).reshape((3,1))
>>> U
matrix([[0],
[0],
[0]])
>>> for k in range(3):
... U[k,0] = randint(-1000,1000)
...
>>> U
matrix([[ 430],
[ -84],
[-316]])
Multiplions 10 fois ce dernier vecteur par $M$
>>> for k in range(10):
... U = M*U
...
>>> U
matrix([[355111006],
[367443780],
[713304158]])
Voici donc une approximation d'un vecteur propre associé à la plus grande valeur propre. Trouvons un plus beau vecteur propre...
>>> U/float(U[0,0])
matrix([[ 1. ],
[ 1.03472935],
[ 2.00867939]])
Donc le vecteur $(1,1,2)^T$ semble bien être un vecteur propre. Vérifions-le, et profitons-en pour trouver ladite plus grande valeur propre...
>>> V = matrix((1,1,2)).reshape((3,1))
>>> M*V
matrix([[-10],
[-10],
[-20]])
>>> M*V==-10*V
matrix([[ True],
[ True],
[ True]], dtype=bool)
On trouve bien -10 pour la plus grande valeur propre (en valeur absolue).
$\left( \begin{array}{ccc} -4 & 4 & -2 \\ 1 & -4 & 5 \\ -4 & 8 & 2 \end{array}\right)$
$\left( \begin{array}{ccc} 15+6i & -16+7i & 5+10i \\ -8+16i & -5-19i & -12+4i \\ 11+23i & -27-11i & -6+18i \end{array}\right)$
On a montré comment obtenir un vecteur propre associé à la valeur propre dominante.
Si la multiplicité de cette valeur propre est plus grande que 1, c'est-à-dire :
alors on peut reproduire l'algorithme plusieurs fois, en partant de vecteurs aléatoires différents : une nouvelle itération a de bonnes chances de conduire à un vecteur propre linéairement indépendant de ceux déjà obtenus.
Pour savoir si les vecteurs obtenus sont bien indépendants, on peut procéder à une élimination de Gauss : si l'on obtient un 0 sur la diagonale du système triangulaire obtenu par la méthode de Gauss, alors la famille n'est pas libre.
On peut montrer que la vitesse de convergence vers la valeur propre dominante $\lambda_1$ est de type exponentiel, avec une raison en général égale à $\frac{\lambda_2}{\lambda_1}$, où $\lambda_2$ est la deuxième plus grande valeur propre (pour le module).
En d'autres termes, le nombre de calculs élémentaires, donc le temps, nécessaire pour obtenir $n$ itérations, est de l'ordre de $\left(\frac{\lambda_2}{\lambda_1}\right)^n$.
Si la matrice n'est pas diagonale, la convergence persiste, mais est bien plus lente (de l'ordre de $\frac{1}{n}$).
Il est possible de trouver toutes les valeurs propres d'une matrice $M$ donnée (toujours telle que ses valeurs propres soient, en module, deux à deux distinctes), en itérant la méthode ci-dessus.
En effet, si :
alors $N = M-\lambda_1 U_1 U_1^T$ possède pour valeurs propres, rangées par module (strictement) décroissant : $(\lambda_2, ..., \lambda_n, 0)$.
On peut donc réappliquer la méthode des puissances itérées à $N$, pour trouver $\lambda_2$, et poursuivre le procédé jusqu'à trouver toutes les valeurs propres.
Notons, cependant, que cette méthode de déflation ne permet pas de trouver la valeur propre nulle.
La méthode précédente permet d'être adaptée au calcul de la plus petite valeur propre $M$ : il suffit de l'appliquer à $M^{-1}$ (ce qui suppose que 0 ne soit pas valeur propre de $M$).
La croissance est rapide, de type exponentielle (de rapport $\frac{|\lambda_1|}{|\lambda_2|}$, où $\lambda_1$ est ici la plus petite valeur propre de $M$, en module, et $\lambda_2$ la deuxième plus petite).
Supposons que l'on connaisse une approximation $\lambda'$ d'une valeur propre $\lambda$ de $M$, et que l'on souhaite améliorer cette connaissance.
On peut montrer que $N=M-\lambda'I_n$ possède une valeur propre égale à $\lambda-\lambda'$, qui est sûrement la plus petite valeur propre de $N$, si $\lambda'$ est une approximation de $\lambda$.
La méthode des puissances itérées inverses permettra donc de trouver $\lambda-\lambda'$, donc $\lambda$ puisque l'on connaît $\lambda'$. De plus, on récupère le vecteur propre $V$, associé à la valeur propre $\lambda-\lambda'$ pour $N$, donc à $\lambda$ pour $M$.
L'intérêt de cette technique est double :
La méthode des puissances itérées inverses est donc, en général, utilisée après une autre méthode donnant seulement une valeur approchée d'une valeur propre : recherche approchée des zéros du polynôme caractéristique, localisation par dichotomie, ou une méthode du texte.
Remarquons, pour finir cette section, que la méthode des puissances itérées inverses appliquée à $N$ suppose d'inverser $N$. Dans la pratique, on évite de le faire :
Dans l'accélération de convergence ci-dessus, on résolvait, à chaque itération, le système $(M-\lambda' I_n) U_{k} = U_{k-1}$, en ayant fait une fois pour toutes la décomposition LU. Dans ce qui précède, $\lambda'$ est l'approximation initiale de la valeur propre recherchée.
La méthode du quotient de Rayleigh consiste à remarquer que $\lambda_k = <M U_k, U_k>$ converge vers $\lambda$ ($<,>$ représente le produit scalaire). On remplace donc la suite de système linéaires $(M-\lambda' I_n) U_{k} = U_{k-1}$ en $(M-\lambda_k I_n) U_{k} = U_{k-1}$
Plus $\lambda_k$ se rapproche de $\lambda$, plus $M-\lambda_k I_n$ risque de ne plus être inversible ($\lambda$ est valeur propre de $M$ si et seulement si $M-\lambda I_n$ n'est pas inversible). La résolution du système $(M-\lambda_k I_n) U_{k} = U_{k-1}$ peut ne plus être possible au-delà d'une certaine itération : il faudra alors stopper le programme.
Pour être exhaustif, finissons ce tour d'horizon des méthodes de détermination des valeurs propres par l'évocation de la méthode dite des réflexions de Householder.
Le principe est ici d'utiliser les matrices de Householder pour ramener la recherche des valeurs propres d'une matrice de taille $n$ à celle d'une matrice de taille $n-1$. Nous n'en dirons pas plus.