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

Answer 2

我认为您的问题是您期望以特定顺序返回特征值，但它们的结果并不总是相同。对它们进行排序，然后您就可以上路了。如果我运行您的代码以使用eig 生成d 和dx，我会得到以下信息：

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