scipy.linalg.eig restituiscono autovalori complessi per la matrice di covarianza?

Question

Gli autovalori di una matrice di covarianza devono essere reali e non negativi perché le matrici di covarianza sono simmetriche e semi positive definite.

Tuttavia, date un'occhiata al seguente esperimento con 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 

Tuttavia, riproducendo l'esempio precedente in Matlab funziona correttamente:

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

lei ha sollevato due questioni:

1. Gli autovalori restituiti da scipy.linalg.eig non sono reali.
2. Alcuni degli autovalori sono negativi.

Entrambi questi problemi sono il risultato di errori introdotti dal troncamento e errori di arrotondamento, che avvengono sempre con algoritmi iterativi utilizzando aritmetica in virgola mobile. Si noti che i risultati Matlab hanno anche prodotto autovalori negativi.

Ora, per un aspetto più interessante del problema: perché il risultato di Matlab è reale, mentre il risultato di SciPy ha alcune componenti complesse?

Matlab's eig rileva se la matrice di input è reale simmetrica o Hermitiana e utilizza la fattorizzazione di Cholesky quando lo è. Vedere la descrizione dell'argomento chol nello eig documentation. Questo non viene fatto automaticamente in SciPy.

Se si desidera utilizzare un algoritmo che sfrutta la struttura di una matrice simmetrica o Hermitiana reale, utilizzare scipy.linalg.eigh. Per l'esempio nella domanda:

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

Questo risultato è lo stesso di Matlab di, se si arrotonda per lo stesso numero di cifre di precisione che Matlab stampato.

Answer 2

cosa si sta verificando è stabilità numerica a causa di limitazioni sulla precisione in virgola mobile.

Nota che:

(1) MATLAB restituito anche valori negativi, ma il formato di stampa è impostato su short e non si vede la precisione completa del doppio memorizzato. Utilizzare format long g per stampare più decimali

(2) Tutte le parti immaginarie restituite da numpy's linalg.eig sono vicine alla precisione della macchina. Quindi dovresti considerarli zero.