L3 MINT M325 Calcul Scientifique II
Université Paris–Sud Centre d’Orsay
AN-EDP Année 2014 − 2015
© J.-B. A. K. <[Link]@[Link]>
Fiche de TP 4 : Méthodes de Lanczos pour calcul de valeurs/vecteurs propres.
On s’intéresse à la détermination des valeurs/vecteurs propres d’une matrice symétrique réelle A, creuse et de grande
taille. Afin de mener à bien la résolution des exercices proposés, on fera usage des utilitaires :
— GSL pour le stockage des matrices et vecteurs.
— Gnuplot pour la représentation graphique des matrices.
Une présentation succincte de ces utilitaires est accessible sur Dokeos dans le dossier dédié au présent cours.
Note 1 (Principe de la méthode de lanczos).
Le principe de la méthode de Lanczos consiste en deux points :
1. Appliquer l’algorithme de Lanczos (qui est une simplification de l’algorithme d’Arnoldi lorsque la
matrice est symétrique) pour construire une base de orthonormée de l’espace de l’espace de Krylov
def
Kk−1 = {r0 , Ar0 , . . . , Ak−1 r0 } ≡ Kk (A, r0 )
Le résultat de cette procédure s’écrit sous la forme :
A Vk = Vk Tk + vk+1 e∗k
où Vk est la matrice dont les colonnes forment la base orthonormée de Kk−1 , Tk est une matrice
tridiagonale symétrique ek est le k-ième vecteur canonique de l’espace Rk , e∗k son adjoint (transposé
dans notre cas réel) et vk+1 le vecteur qui, normalisé, serait le prochain élément de la base de Kk .
2. Calculer les valeurs et vecteurs propres de la matrice tridiagonale Tk (soit analytiquement soit par
une méthode de dichotomie dite de Givens). Et déduire les approximations des vecteurs et valeurs
propres de A.
On constate que pour cette méthode :
— la matrice A n’intervient que sous forme de produit matrice-vecteur ;
— Vk peut être stockée comme une matrice de taille n × k ou comme un tableaux de k vecteurs de taille
n chacun.
— Tk étant tridiagonale et symétrique, elle peut être stockée en utilisant deux vecteurs ; un de taille k
pour la diagonale principale et un autre de taille k − 1 pour la première sur-diagonale.
Thème - 1 Méthode de bissection (ou dichotomie) de Givens
Soit A ∈ Mn (R) une matrice tridiagonale symétrique.
b0 c0 0
c0 . . . ..
.
A= ,
.. ..
. . cn−2
0 cn−2 bn−1
où ci 6= 0, ∀0 ≤ i ≤ n − 1.
1
Q-1 : On désigne par Ai la sous-matrice principale d’ordre i de A et on pose pi (λ) = det(Ai −λIi ), son polynôme
caractéristique.
Q-1-1 : En développant le déterminant det(Ai − λIi ) par rapport à la dernière colonne, montrer que la suite pi
vérifie la relation de récurrence
p0 (λ) = 1, p1 (λ) = b0 − λ,
pi (λ) = (bi−1 − λ)pi−1 (λ) − c2i−2 pi−2 (λ), ∀i ≥ 2.
Q-1-2 : Écrire une fonction C de prototype
double pc(const gsl vector* b, const gsl vector* c, int i, double lambda)
définissant la i−ème fonction de cette suite, où b et c sont respectivement des vecteurs de taille n et n − 1 représentant
la diagonale principale et la première sur-diagonale de la matrice tridiagonale symétrique A.
Q-1-3 : Écrire une fonction C de prototype
int Nracines(const gsl vector* b, const gsl vector* c, int i, double mu)
qui calcule pour une matrice tridiagonale symétrique A (stockée dans les vecteurs b et c), le nombre de racines de pi
qui sont strictement inférieures à mu, comme étant le nombre de changement de signe entre éléments consécutifs de
l’ensemble {1, sgnp1 (µ), . . . sgnpi (µ)}, où on a posé
signe de pi (µ) si pi (µ) 6= 0,
sgnpi (µ) =
signe de pi−1 (µ) si pi (µ) = 0.
Q-2 : On suppose que la matrice A a n valeurs propres distinctes λ0 < . . . < λn−1 .
Q-2-1 : Soit 0 ≤ i ≤ n − 1 un entier. On suppose que la valeur propre λi se trouve dans l’intervalle [r0 , s0 ].
Montrer que :
Nracines(A, n, r0 +s
2 ) ≥ i + 1 =⇒ λi ∈ [r0 ,
0 s0 +r0
2 [
r0 +s0 r0 +s0
Nracines(A, n, 2 ) < i + 1 =⇒ λi ∈ [ 2 , s0 ].
En déduire qu’on peut déterminer la valeur propre λi par dichotomie avec une précision ε donnée.
Q-2-2 : Ecrire une fonction C de prototype double bissectionGivens(const gsl vector* b,
const gsl vector* c, int i, double r0, double s0, double tol) qui détermine la i−ème va-
leur propre de la matrice tridiagonale A (stockée sous forme bande dans les vecteurs b et c), avec une précision tol, où
[a0,b0] est un intervalle de départ contenant ladite valeur propre (on pourra prendre r0 = −kAk∞ ,s0 = kAk∞ ).
2
Thème - 2 Méthode de Lanczos
Soit A ∈ Mn (R) une matrice symétrique.
Q-3 : Programmer l’algorithme de Lanczos à travers une fonction C
Prototype de fonction pour l’algorithme de Lanczos
void lanczos(const MatriceCSR* A, gsl_vector* r0, int k, gsl_matrix* V, gsl_matrix* T, int* k0)
qui prend comme arguments
• la matrice symétrique A,
• un vecteur r0 et un entier k définissant l’espace de Krylov Kk = {r0 , Ar0 , . . . , Ak r0 },
et fournit en retour
• un entier k0 égal à k ou coı̈ncidant avec la dimension critique de Krylov si k est supérieur à la dimension critique
de Krylov.
• une matrice V de taille n × k0 dont les colonnes sont les k0 premiers vecteurs de la base orthonormée de l’espace
de Krylov Kk0 ,
• une matrice T tridiagonale de taille k0 × k0 telle que V ∗ AV = T .
L’algorithme en pseudo code est de la forme
Algorithme de Lanczos
----------------------------------------------------------------------------------------------------
ALGORITHME DE LANCZOS | COMMENTAIRES
---------------------------------------------------------------------------------------------------
Donnees: A, r0, k |
Resultat V, T, k0 |
-----------------------------------------------|---------------------------------------------------
Initialization : |
|
w = r0/|| r_0 || |
T mis a 0 |
V(:,0) = w | GSl dispose d’une fonction g s l m a t r i x s e t c o l
for j = 0 à k-1 |
w = A * V(:,j); | GSl dispose d’une fonction g s l m a t r i x g e t c o l
for i = max(0,j-1) à j |
T(i,j) = (w, V(:,i)); |
w = w - T(i,j) * V(:,i); |
end |
normj = ||w||; T(j+1,j) = normj |
|
if(normj != 0 ) | Il vaut mieux faire: normj > epsilon
V(:,j+1) = w/normj | GSl dispose d’un fonction g s l m a t r i x s e t c o l
else |
k0 = j STOP. |
return |
end |
end |
k0 = k; |
----------------------------------------------------------------------------------------------------
Pour la programmation en C , T arrive avec une taille kmax x kmax, et V avec une taille n x kmax
En sortie de la fonction, comme on dispose de k0, on peut redimensionner V et T. ( kmax = k+1 )
Q-4 :
Q-4-1 : Générer aléatoirement une matrice symétrique définie positive A de taille n × n, et déterminer sa plus
grande valeur propre λA à l’aide de la fonction (des précédents Tps) valeurs propres extremes.
Q-4-2 : Faire varier k = 2, . . . n, et comparer à chaque fois la plus grande valeur propre de T obtenue par la
méthode de bissection de l’exercice précédent à λA .
(On pourra prendre pour r0 le premier vecteur colonne de la matrice A).
Q-4-3 : A-t-on toujours AV = V T ? Comment évolue la valeur absolue de la différence des plus grandes valeurs
propres de A et de T en fonction de la norme de Frobenius de AV − V T ? En fonction de la dernière composante du
vecteur propre de T associé à sa plus grande valeur propre ?