Bonjour,
Je cherche depuis hier une solution à mon problème de calcul de valeurs propres:
J'ai une matrice plutôt imposante de taille ~(5000, 5000) et il faut que je calcule ses valeurs propres et ses vecteurs propres. Par construction, la matrice est symétrique donc les calculs sont optimisés! Néanmoins, les résultats donnés par le code ci-dessous ne me suffisent pas, il me semblent trop 'approchés' de la valeur vraie:
sgemm étant la fonction de multiplication matricielle de BLAS 3. Elle ne pose aucun problème, le programme réalise en amont des opérations bien plus grandes avec cette fonction sans qu'aucun problème ne soit survenu (une quinzaine de multiplications matricielles très grandes: (5278, 49013)x(49013, 5278)).
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41 import scipy.linalg import numpy from scipy.linalg.blas import sgemm # Création de la matrice à décomposer, on va l'appeler G. E1, V1 = scipy.linalg.eig(G) E1s = numpy.sort(E1) V1s = numpy.sort(V1) E2, V2 = scipy.linalg.eigh(G) E2s = numpy.sort(E2) V2s = numpy.sort(V2) E3 = scipy.linalg.eigvals(G) E3s = numpy.sort(E3) E4 = scipy.linalg.eigvalsh(G) E4s = numpy.sort(E4) # On vérifie le résultat: print( "Écarts entre les valeurs propres:", "\n\t||E2-E1||:", scipy.linalg.norm(E2s-E1s), "\n\t||E3-E1||:", scipy.linalg.norm(E3s-E1s), "\n\t||E4-E1||:", scipy.linalg.norm(E4s-E1s), "\n\t||E3-E2||:", scipy.linalg.norm(E3s-E2s), "\n\t||E4-E2||:", scipy.linalg.norm(E4s-E2s), "\n\t||E4-E3||:", scipy.linalg.norm(E4s-E3s)) print( "Écarts entre les vecteur propres:", "\n\t||V2-V1||:", fNorm(V2s-V1s)) def fNorm(A): # Retourne la racine carrée de la somme de tous les éléments au carré. Une norme Euclidienne, mais sur des matrices :) return np.sum(A**2)**0.5 def verify(Initial, Values, Vectors): # fNorm(A) -> retourne la norme de Frobenius de A # sgemm(alpha, A, B) = alpha*(AxB) = produit matriciel de A par B, multiplié par un coefficient alpha return fNorm(sgemm(alpha=1.0, a=Initial, b=Vectors) - sgemm(alpha=1.0, a=Vectors, b=numpy.diag(Values))) print( "Écarts entre les résultats et ce qui est attendu:", "\n\t||G*V1-V1*D(E1)||:", verify(G, E1, V1), "\n\t||G*V1-V1*D(E3)||:", verify(G, E3, V1), "\n\t||G*V2-V2*D(E2)||:", verify(G, E2, V2), "\n\t||G*V2-V2*D(E4)||:", verify(G, E4, V2))
Et voici la sortie:
Mettons à part le fait que eig et eigh ne donnent pas du tout les mêmes valeurs, c'est surtout les dernières lignes qui m’intéressent: l'erreur calculée est au minimum de 2.63! En comparaison, le même code en Matlab, avec la même matrice me donne une erreur de l'ordre de 1e-7!
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14 Écarts entre les valeurs propres: ||E2-E1||: 87.20174407958984 ||E3-E1||: 7.196725368499756 ||E4-E1||: 101.92095947265625 ||E3-E2||: 87.2437515258789 ||E4-E2||: 19.32546615600586 ||E4-E3||: 101.96918487548828 Écarts entre les vecteur propres: ||V2-V1||: 0.0611834875335 Écarts entre les résultats et ce qui est attendu: ||G*V1-V1*D(E1)||: 43.2232758902 ||G*V1-V1*D(E3)||: 183.425095642 ||G*V2-V2*D(E2)||: 2.72478850655 ||G*V2-V2*D(E4)||: 2.64568433509
J'ai aussi utilisé les fonctions LAPACK via Scipy, et il se trouve que si j'utilise DSYEVR (double précision), l'erreur obtenue est divisée par 2, mais ça reste largement en dessous de Matlab
Comment est-ce que je pourrais palier à ce problème de précision?
Partager