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

特征值分解满足 A V = V Lambda，这是所有可以保证的——例如特征值的顺序不是。

回答你问题的第二部分:

现代编译器/线性代数库产生/包含做不同事情的代码
取决于数据是否在内存中对齐(例如)16 字节边界。这会影响计算中的舍入误差，因为浮点运算以不同的顺序进行。如果算法(此处为 LAPACK/xGEEV)不能保证这方面的数值稳定性，则舍入误差的微小变化可能会影响特征值的排序。

(如果你的代码对这样的事情很敏感，这是不正确的!例如在不同的平台或不同的库版本上运行会导致类似的问题。)

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

见 http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf更多