numpy/scipy eigendecompositions

Question

के लिए अनियमित परिणाम मुझे लगता है कि scipy.linalg.eig कभी-कभी असंगत परिणाम देता है। लेकिन हर बार नहीं।

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

जब मैं बड़ी रेखीय बीजगणित जादूगर नहीं हूँ, मुझे समझ नहीं है कि eigendecomposition स्वाभाविक अजीब राउंडिंग त्रुटियों के अधीन है, लेकिन मुझे समझ नहीं आता क्यों गणना दोहरा भिन्न मान जाएगी। लेकिन मेरे परिणाम और पुनरुत्पादन भिन्न हैं।

समस्या की प्रकृति वास्तव में क्या है - ठीक है, कभी-कभी परिणाम स्वीकार्य रूप से भिन्न होते हैं, और कभी-कभी वे नहीं होते हैं। यहां कुछ उदाहरण दिए गए हैं:

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

~ 3e-13 के ऊपर का अंतर एक बहुत बड़ा सौदा नहीं लगता है। इसके बजाए, असली समस्या (कम से कम मेरे वर्तमान प्रोजेक्ट के लिए) यह है कि कुछ ईजीनवे उचित संकेत पर सहमत नहीं लग सकते हैं।

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

MATLAB में इसी तरह की कोड में यह समस्या नहीं प्रतीत होती है।

Answer 1

eigenvalue decompositions ए वी संतुष्ट = वी लैम्ब्डा, जो सभी क्या गारंटी है है --- उदाहरण के लिए eigenvalues ​​का आदेश नहीं है। अपने प्रश्न के दूसरे भाग के लिए

उत्तर:

आधुनिक compilers/रेखीय बीजगणित पुस्तकालयों का उत्पादन/कि अलग अलग बातें डेटा (जैसे) में 16-बाइट सीमाओं पर स्मृति में गठबंधन है जो इस पर निर्भर करता है कोड होता है। यह कंप्यूटेशंस में गोल करने वाली त्रुटि को प्रभावित करता है, क्योंकि फ़्लोटिंग पॉइंट ऑपरेशंस एक अलग क्रम में किया जाता है। राउंडिंग त्रुटि में छोटे बदलाव तब इग्निग्लूम (यहां, LAPACK/xGEEV) इस संबंध में संख्यात्मक स्थिरता की गारंटी नहीं देते हैं, जैसे eigenvalues ​​को ऑर्डर करने जैसी चीजों को प्रभावित कर सकते हैं।

(यदि आपका कोड इस तरह बातें के प्रति संवेदनशील है, यह सही नहीं है एक अलग मंच या अलग पुस्तकालय संस्करण पर चल रहा है जैसे यह एक ऐसी ही समस्या पैदा होता!।)

परिणामों को आम तौर अर्ध नियतात्मक हैं - - उदाहरण के लिए आपको 2 संभावित परिणामों में से एक मिलता है, भले ही सरणी स्मृति में गठबंधन हो या नहीं। यदि आप संरेखण के बारे में उत्सुक हैं, तो A.__array_interface__['data'][0] % 16 देखें।

के लिए और अधिक

Answer 2

मुझे लगता है कि आपकी समस्या यह है कि आप किसी विशेष क्रम में ईगेंवल्यू वापस लौटने की उम्मीद कर रहे हैं, और वे हमेशा बाहर नहीं आते हैं। उन्हें क्रमबद्ध करें, और आप अपने रास्ते पर होंगे। मैं निम्नलिखित मिल अगर मैं अपने कोड और ddxeig साथ उत्पन्न करने के लिए चलाएँ:

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

... लेकिन

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