2016-04-24 42 views
5

Pozitif semidefinite matrisler için negatif özdeğerler döndüren scipy'nin eigh işleviyle ilgili bazı sorunlar yaşıyorum. Aşağıda bir MWE.scipy eigh pozitif semidefinite matrisi için negatif özdeğerler verir

hess_R işlevi, pozitif bir semidefinite matrisini döndürür (bu, negatif olmayan girdilerle birlikte bir sıra matrisinin ve çapraz matrisin toplamıdır). Basılı

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 

özdeğerler ben hess_R dönüş açıklamada np.float32 ile np.float64 değiştirin yerine

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

olsun

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

, bu yüzden bu tahmin ediyorum çeşit hassas konu.

Bunu düzeltmenin bir yolu var mı? Teknik olarak benim için kullanmam gerekmiyor, ama bence bu benim diğer hatalarla (bu matrislerin kareköklerini alarak, NaN almak, vb.) Altta yatan problem olduğunu düşünüyorum.

+0

"LA.eigh" yerine "LA.eig" yi kullanırsam farklı özdeğerler alırım: '[5.15260753e + 18 + 0.j 3.22785571e + 01 + 0.j 4.62855397e + 16 + 0.j ] – Peaceful

+2

IMHO, Hess_R işleviniz gerçek bir Hessian Matrix döndürmez. yani 'geri çekil' davasında yanlış sonuç döndürüyor. –

+0

@ B.M. Ne demek istediğini açıklayabilir misin? Bunun yerine dönen işlev nedir? – angryavian

cevap

0

Sanırım sorun siz kayan nokta hassasiyeti sınırları. Doğrusal cebir sonuçları için iyi bir kural, float32 için 10^8'de bir parça ve float 64 için 10^16'da yaklaşık bir parça için iyi olmalarıdır. En küçükten en büyüğüne kadar olan orana benziyor. Burada özdeğer 10^-16'dan azdır. Bu nedenle, döndürülen değer gerçekten güvenilir olamaz ve kullandığınız özdeğer uygulamasının ayrıntılarına bağlı olacaktır. Örneğin, burada sahip olmanız gereken dört farklı çözücü vardır; Onların sonuçlarına bir göz atın:

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

Çıktı: hepsi büyük iki özdeğerler üzerinde anlaşmaya

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

Bildirim, ama bu uygulamaya uygulamadan en küçük özdeğer değeri değişir. Yine de, dört durumda da, giriş matrisi 64-bit hassasiyete kadar yeniden yapılandırılabilir: bu, algoritmaların, kendilerine sunulan hassasiyete kadar beklendiği şekilde çalıştığı anlamına gelir.