scipy.linalg.eig zwracają złożone wartości własne macierzy kowariancji?

Question

Wartości własne macierzy kowariancji powinny być rzeczywiste i nieujemne, ponieważ macierze kowariancji są symetryczne i półznacznie określone.

jednak przyjrzeć się następującym doświadczeniu z scipy:

>>> a=np.random.random(5) 
>>> b=np.random.random(5) 
>>> ab = np.vstack((a,b)).T 
>>> C=np.cov(ab) 
>>> eig(C) 
7.90174997e-01 +0.00000000e+00j, 
2.38344473e-17 +6.15983679e-17j, 
2.38344473e-17 -6.15983679e-17j, 
-1.76100435e-17 +0.00000000e+00j, 
5.42658040e-33 +0.00000000e+00j 

Jednak, odtwarzając powyższy przykład w Matlab działa poprawnie:

a = [0.6271, 0.4314, 0.3453, 0.8073, 0.9739] 
b = [0.1924, 0.3680, 0.0568, 0.1831, 0.0176] 
C=cov([a;b]) 
eig(C) 
-0.0000 
-0.0000 
0.0000 
0.0000 
0.7902 

Answer 1

pan podniósł dwie kwestie:

1. Do wartości własne zwracane przez scipy.linalg.eig nie są prawdziwe.
2. Niektóre wartości własne są ujemne.

Oba te problemy są wynikiem błędów wprowadzonych przez błędy skracania i zaokrąglania, co zawsze ma miejsce w przypadku algorytmów iteracyjnych wykorzystujących arytmetykę zmiennoprzecinkową. Zauważ, że wyniki Matlaba również generują ujemne wartości własne.

Teraz, dla bardziej interesującego aspektu problemu: dlaczego wynik Matlab jest realny, podczas gdy wynik SciPy ma pewne złożone elementy?

Matlab's eig wykrywa, czy matryca wejściowa jest prawdziwie symetryczna lub hermetyzowana, i używa Cholesky factorization, gdy jest. Zobacz opis argumentu chol w eig documentation. Nie odbywa się to automatycznie w SciPy.

Jeśli chcesz użyć algorytmu wykorzystującego strukturę rzeczywistej macierzy symetrycznej lub hermitowskiej, użyj scipy.linalg.eigh. Na przykład w pytaniu:

>>> eigh(C, eigvals_only=True) 
array([ -3.73825923e-17, -1.60154836e-17, 8.11704449e-19, 
     3.65055777e-17, 7.90175615e-01]) 

Ten wynik jest taki sam jak na Matlab, jeśli zaokrąglić do tej samej liczby cyfr precyzji że Matlab drukowanych.

Answer 2

czego doświadczamy jest niestabilność numeryczna ze względu na ograniczenia precyzji zmiennoprzecinkowej.

zauważyć, że:

(1) MATLAB powrócił także wartości ujemne, ale druk jest ustawiony na short i nie zobaczyć pełną precyzję podwójne przechowywane w pamięci. Użyj format long g do drukowania większej liczby miejsc po przecinku.

(2) Wszystkie wyimaginowane części zwrócone przez linalg.eig z numpy są zbliżone do precyzji maszyny. Dlatego powinieneś rozważyć je jako zero.