scipy eigh fornisce autovalori negativi per matrice semidefinita positiva

Question

Sto riscontrando alcuni problemi con la funzione eigh di scipy che restituisce autovalori negativi per matrici semidefinite positive. Di seguito è riportato un MWE.

La funzione hess_R restituisce una matrice semidefinita positiva (è la somma di una matrice di rango uno e di una matrice diagonale, entrambe con voci non negative).

import numpy as np 
from scipy import linalg as LA 

def hess_R(x): 
    d = len(x) 
    H = np.ones(d*d).reshape(d,d)/(1 - np.sum(x))**2 
    H = H + np.diag(1/(x**2)) 
    return H.astype(np.float64) 

x = np.array([ 9.98510710e-02 , 9.00148922e-01 , 4.41547488e-10]) 
H = hess_R(x) 
w,v = LA.eigh(H) 
print w 

Gli autovalori stampati sono

[ -6.74055241e-271 4.62855397e+016 5.15260753e+018] 

Se sostituisco np.float64 con np.float32 nella dichiarazione ritorno di hess_R ottengo

[ -5.42905303e+10 4.62854925e+16 5.15260506e+18] 

invece, in modo da sto indovinando questo è una sorta di problema di precisione.

C'è un modo per risolvere questo problema? Tecnicamente non ho bisogno di usare Eigh, ma penso che questo sia il problema di fondo con i miei altri errori (prendendo radici quadrate di queste matrici, ottenendo NaNs, ecc.)

Answer 1

Penso che il problema è che hai colpito il limiti di precisione in virgola mobile. Una buona regola per i risultati di algebra lineare è che sono validi per circa una parte in 10^8 per float32 e circa una parte in 10^16 per float 64. Sembra che il rapporto tra il più piccolo e il più grande l'autovalore qui è inferiore a 10^-16. Per questo motivo, il valore restituito non può essere considerato attendibile e dipenderà dai dettagli dell'implementazione autovalore che si utilizza.

Ad esempio, ecco quattro diversi risolutori che dovrebbero essere disponibili; dare un'occhiata al loro risultati:

# using the 64-bit version 
for impl in [np.linalg.eig, np.linalg.eigh, LA.eig, LA.eigh]: 
    w, v = impl(H) 
    print(np.sort(w)) 
    reconstructed = np.dot(v * w, v.conj().T) 
    print("Allclose:", np.allclose(reconstructed, H), '\n') 

uscita:

[ -3.01441754e+02 4.62855397e+16 5.15260753e+18] 
Allclose: True 

[ 3.66099625e+02 4.62855397e+16 5.15260753e+18] 
Allclose: True 

[ -3.01441754e+02+0.j 4.62855397e+16+0.j 5.15260753e+18+0.j] 
Allclose: True 

[ 3.83999999e+02 4.62855397e+16 5.15260753e+18] 
Allclose: True 

Avviso sono tutti d'accordo sulle grandi due autovalori, ma che il valore dei piccoli cambiamenti agli autovalori da implementazione a implementazione. Tuttavia, in tutti e quattro i casi la matrice di input può essere ricostruita con precisione fino a 64 bit: questo significa che gli algoritmi funzionano come previsto fino alla precisione a loro disposizione.