numpy/scipy特征分解的不稳定结果

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

虽然我不是最好的线性代数向导,但我确实理解特征分解固有地受到奇怪的舍入误差的影响,但我不明白为什么重复计算会导致不同的值.但我的结果和可重复性是不同的.

问题的本质究竟是什么 - 好吧,有时结果是可以接受的,有时它们不是.这里有些例子:

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

特征值分解满足AV = V Lambda,这是所有保证的 - 例如特征值的顺序不是.

回答你问题的第二部分:

现代编译器/线性代数库产生/包含代码,这些代码根据数据是否在(例如)16字节边界的内存中对齐而执行不同的操作.这会影响计算中的舍入误差,因为浮点运算以不同的顺序完成.如果算法(这里,LAPACK/xGEEV)在这方面不保证数值稳定性,则舍入误差的微小变化可以影响诸如特征值的排序之类的事物.

(如果您的代码对此类内容敏感,则不正确!在不同平台或不同库版本上运行会导致类似问题.)

结果通常是准确定性的 - 例如,你得到2个可能的结果之一,这取决于数组是否恰好在内存中对齐.如果您对对齐感到好奇,请检查A.__array_interface__['data'][0] % 16.

有关更多信息,请参见http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf