2017-09-28 3 views
2

Gibt es eine Möglichkeit, die Genauigkeit der Ausgabe von numpy.linalg.eig() und scipy.linalg.eig() zu verbessern?Rundungsfehler in numpy.linalg.eig() und scipy.linalg.eig()

Ich diagonale eine nicht-symmetrische Matrix, aber ich erwarte aus physikalischen Gründen ein reales Spektrum von Paaren von positiven und negativen Eigenwerten. Tatsächlich kommen die Eigenwerte paarweise, und ich habe durch eine unabhängige analytische Berechnung verifiziert, dass zwei der Paare richtig sind. Das problematische Paar ist dasjenige mit Eigenwerten nahe Null, die kleine Imaginärteile zu haben scheinen. Ich erwarte, dass dieses Paar bei Null entartet ist, so dass die Imaginärteile höchstens maschinell präzise sein können, aber sie sind viel größer. Ich denke, dass dies zu einem kleinen Fehler in den Eigenvektoren führt, der sich jedoch in nachfolgenden Manipulationen ausbreiten kann.

Das folgende Beispiel zeigt, dass es fiktive Imaginärteil-Reste gibt, indem die Gültigkeit der Transformation überprüft wird.

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

die

3.93435308362e-09 

eine Anzahl von der Ordnung des fiktiven imaginären Teils des Nulleigenwertes gibt.

+0

erwähnt Wie in den Antworten, Invertierung a:

Aus diesem Satz, in den schlechteren Umständen ein Fehler bei Maschinengenauigkeit im Eingangsbereich zu einem Fehler in der Größenordnung von 1e-8 in den Eigenwerten seit ausbreiten könnte Matrix numerisch führt normalerweise zu schlechten Ergebnissen. Es gibt fast immer einen besseren Weg, um Ihr Problem zu lösen, ohne eine Matrix zu invertieren. –

Antwort

2

Sie verlieren möglicherweise etwas Genauigkeit, indem Sie die Inverse beim Berechnen des obigen Fehlers verwenden. Wenn Sie stattdessen berechnen:

# 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 

der Fehler ist viel kleiner.

Ihr Problem kann auch unter schlechter Konditionierung leiden, was bedeutet, dass es eine sehr hohe numerische Empfindlichkeit für Rundungen und andere Fehler geben wird. Die Bauer-Fike Theorem gibt eine Obergrenze für die Fehlerempfindlichkeit des Eigenwertproblems.

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