5

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

私は最高の線形代数の魔法使いではありませんが、固有分解が本質的に奇妙な丸め誤差の影響を受けることは理解していますが、計算を繰り返すと異なる値になる理由がわかりません。しかし、私の結果と再現性はさまざまです。

問題の性質は正確には何ですか - まあ、結果が許容できるほど異なる場合もあれば、そうでない場合もあります。ここではいくつかの例を示します。

>>> 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 の同様のコードには、この問題はないようです。

4

2 に答える 2

5

固有値分解は、AV = V ラムダを満たします。これはすべて保証されているものです。たとえば、固有値の順序はそうではありません。

質問の 2 番目の部分への回答:

最新のコンパイラ/線形代数ライブラリは、データが (たとえば) 16 バイト境界でメモリ内で整列されているかどうかに応じて、さまざまな処理を行うコードを生成/含んでいます。これは、浮動小数点演算が異なる順序で行われるため、計算の丸め誤差に影響します。丸め誤差の小さな変化は、アルゴリズム (ここでは LAPACK/xGEEV) がこの点で数値の安定性を保証しない場合、固有値の順序付けなどに影響を与える可能性があります。

(あなたのコードがこのようなことに敏感なら、それは正しくありません!例えば、別のプラットフォームや別のライブラリ バージョンでコードを実行すると、同様の問題が発生します。)

通常、結果は準決定論的です --- たとえば、配列がたまたまメモリ内でアラインされているかどうかに応じて、2 つの可能な結果のうちの 1 つが得られます。アライメントに興味がある場合は、 を確認してくださいA.__array_interface__['data'][0] % 16。

詳細については、 http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdfを参照してください。

于 2013-02-04T23:26:44.733 に答える
2

あなたの問題は、固有値が特定の順序で返されることを期待していることであり、それらが常に同じになるとは限らないと思います。それらを並べ替えれば、あなたはあなたの道にいるでしょう。コードを実行して生成するdとdx、eig次のようになります。

>>> 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)
于 2013-02-04T21:00:39.543 に答える