risultati erratici per eigendecompositions numpy/scipy

Question

Sto trovando che scipy.linalg.eig a volte dà risultati incoerenti. Ma non ogni volta.

>>> import numpy as np 
>>> import scipy.linalg as lin 
>>> modmat=np.random.random((150,150)) 
>>> modmat=modmat+modmat.T # the data i am interested in is described by real symmetric matrices 
>>> d,v=lin.eig(modmat) 
>>> dx=d.copy() 
>>> vx=v.copy() 
>>> d,v=lin.eig(modmat) 
>>> np.all(d==dx) 
False 
>>> np.all(v==vx) 
False 
>>> e,w=lin.eigh(modmat) 
>>> ex=e.copy() 
>>> wx=w.copy() 
>>> e,w=lin.eigh(modmat) 
>>> np.all(e==ex) 
True 
>>> e,w=lin.eigh(modmat) 
>>> np.all(e==ex) 
False 

Mentre io non sono il più grande mago di algebra lineare, capisco che l'eigendecomposition è intrinsecamente soggetto a errori di arrotondamento strani, ma non capisco il motivo per ripetere il calcolo si tradurrebbe in un valore diverso. Ma i miei risultati e riproducibilità sono variabili.

Qual è esattamente la natura del problema - beh, a volte i risultati sono accettabilmente diversi e, a volte, non lo sono. Ecco alcuni esempi:

>>> d[1] 
(9.8986888573772465+0j) 
>>> dx[1] 
(9.8986888573772092+0j) 

La differenza sopra di ~ 3e-13 non sembra un grande affare enorme. Invece, il vero problema (almeno per il mio progetto attuale) è che alcuni degli autovalori non sembrano concordare sul segno corretto.

>>> np.all(np.sign(d)==np.sign(dx)) 
False 
>>> np.nonzero(np.sign(d)!=np.sign(dx)) 
(array([ 38, 39, 40, 41, 42, 45, 46, 47, 79, 80, 81, 82, 83, 
    84, 109, 112]),) 
>>> d[38] 
(-6.4011617320002525+0j) 
>>> dx[38] 
(6.1888785138080209+0j) 

Codice simile a MATLAB non sembra presentare questo problema.

Answer 1

Le scomposizioni autovalori soddisfano AV = V Lambda, che è tutto ciò che è garantito --- per esempio l'ordine degli autovalori non lo è.

risposta alla seconda parte della tua domanda:

compilatori moderni/librerie di algebra lineare producono/contengono codice che fa cose diverse a seconda che i dati vengono allineati in memoria (per esempio) i limiti di 16 byte. Ciò influenza l'errore di arrotondamento nei calcoli, poiché le operazioni in virgola mobile vengono eseguite in un ordine diverso. Piccole variazioni nell'errore di arrotondamento possono quindi influire su cose come l'ordinamento degli autovalori se l'algoritmo (qui, LAPACK/xGEEV) non garantisce la stabilità numerica a questo riguardo.

(Se il codice è sensibile a cose come questa, che non è corretto in esecuzione ad esempio su una piattaforma diversa o una versione libreria diversa porterebbe a un problema simile!).

I risultati di solito sono quasi-deterministico - - Ad esempio, si ottiene uno dei 2 possibili risultati, a seconda che l'array sia allineato in memoria o meno. Se sei curioso di allineare, controlla A.__array_interface__['data'][0] % 16.

Vedi http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf per più

Answer 2

Penso che il tuo problema è che ti aspetti che gli autovalori vengano restituiti in un ordine particolare, e non sempre escono allo stesso modo. Ordinali, e sarai sulla buona strada. Se corro il codice per generare d e dx con eig ottengo il seguente:

>>> np.max(d - dx) 
(19.275224236664116+0j) 

Ma ...

>>> d_i = np.argsort(d) 
>>> dx_i = np.argsort(dx) 
>>> np.max(d[d_i] - dx[dx_i]) 
(1.1368683772161603e-13+0j)