Discussion :

[Invité]

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:

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))

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)).

Et voici la sortie:

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

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!
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?

[__dardanos__]

Salut,
Les vecteurs propres retournés par les fonctions eig et eigh ne sont peut-être pas classés dans le même ordre. Ça me semble risqué de directement faire E2 - E1.
On ne pas non plus savoir si les fonctions sgemm et fNorm sont correctes (d'où viennent-elles ?).

[Invité]

J'ai édité mon premier message afin de faciliter la lecture d'éventuels nouveaux arrivants Tu peux trouver toutes les infos que tu cherches dans ce message