2017-09-28 3 views
2

Y at-il un moyen d'améliorer la précision de la sortie de numpy.linalg.eig() et scipy.linalg.eig()?erreurs d'arrondi dans numpy.linalg.eig() et scipy.linalg.eig()

Je diagonise une matrice non-symétrique, mais je m'attends à ce que, physiquement, j'obtienne un spectre réel de paires de valeurs propres positives et négatives. En effet, les valeurs propres viennent par paires, et j'ai vérifié par un calcul analytique indépendant que deux des paires sont correctes. La paire problématique est celle avec des valeurs propres proches de zéro, qui semblent avoir de petites parties imaginaires. Je m'attends à ce que cette paire soit dégénérée à zéro, donc les parties imaginaires peuvent tout au plus être à la précision de la machine, mais elles sont beaucoup plus grandes. Je pense que cela conduit à une petite erreur dans les vecteurs propres, qui peuvent cependant se propager dans des manipulations ultérieures.

L'exemple ci-dessous montre qu'il existe des restes fictifs imaginaires, en vérifiant la validité de la transformation.

import numpy as np 
import scipy.linalg as sla 

H = np.array(
    [[ 11.52, -1., -1.,  9.52, 0.,  0. ], 
     [ -1., 11.52, -1.,  0.,  9.52, 0., ], 
     [ -1., -1., 11.52, 0.,  0.,  9.52,], 
     [ -9.52, 0.,  0., -11.52, 1.,  1., ], 
     [ 0., -9.52, 0.,  1., -11.52, 1., ], 
     [ 0.,  0., -9.52, 1.,  1., -11.52 ]], 
    dtype=np.float64 
      ) 

#print(H) 
E,V = np.linalg.eig(H) 
#E,V = sla.eig(H) 
H2=reduce(np.dot,[V,np.diag(E),np.linalg.inv(V)]) 
#print(H2) 
print(np.linalg.norm(H-H2)) 

qui donne

3.93435308362e-09 

un certain nombre de l'ordre de la partie imaginaire fictive de zéro des valeurs propres.

+0

Comme mentionné dans les réponses, inversion d'un matrice numériquement conduira généralement à de mauvais résultats. Il y a presque toujours une meilleure façon de résoudre votre problème sans inverser une matrice. –

Répondre

2

Vous pouvez perdre de la précision en prenant l'inverse en calculant l'erreur ci-dessus. Si à la place vous calculez:

# H = V.E.inv(V) <=> H.V = V.E 
print(np.linalg.norm(H.dot(V)-V.dot(np.diag(E)))) 
# error: 2.81034671113e-14 

l'erreur est beaucoup plus petite.

Votre problème peut également souffrir d'un mauvais conditionnement, ce qui signifie qu'il y aura une très grande sensibilité numérique aux arrondis et autres erreurs. Le Bauer-Fike Theorem donne une limite supérieure sur la sensibilité aux erreurs du problème de valeur propre.

De ce théorème, dans les circonstances pire une erreur à la précision de la machine dans le domaine d'entrée pourrait se propager à une erreur de l'ordre de 1e-8 dans les valeurs propres depuis:

machine_precison = np.finfo(np.float64).eps 
print(np.linalg.cond(V)*(machine_precison)) 
# 4.54517272701e-08