Numpy/Scipy eigh 函数的结果不一致
Inconsistent Results for Numpy/Scipy eigh function
我从 scipy.linalg.eigh 函数中得到不一致的结果（如果我尝试使用 numpy.linalg.eigh 也是如此）。即使我使用相同的输入，我得到的值也是不同的。为了以防万一，我尝试将随机种子设置为一个常量，但还是不行。
这是我的测试代码：
print('#'*50)
import numpy as np
import scipy
from scipy.linalg import eigh
# Make a random symmetric matrix
m_size = 1500
A = scipy.random.rand(m_size, m_size)
L = np.triu(A) + np.triu(A,1).T
# Verify that matrix is symmetric
print('L symmetric: %s' % np.array_equal(L,L.T))
print('')
# Just to make sure that changes to one L won't effect changes to another
# and that we are starting with the exact same input.
import copy
L1 = copy.deepcopy(L)
L2 = copy.deepcopy(L)
L3 = copy.deepcopy(L)
L4 = copy.deepcopy(L)
L5 = copy.deepcopy(L)
# Sanity check that inputs are the same
print('L == L1 (before): %s' % str(np.array_equal(L,L1)))
print('L == L2 (before): %s' % str(np.array_equal(L,L2)))
print('L == L3 (before): %s' % str(np.array_equal(L,L3)))
print('L == L4 (before): %s' % str(np.array_equal(L,L4)))
print('L == L5 (before): %s' % str(np.array_equal(L,L5)))
print('L1 == L2 (before): %s' % str(np.array_equal(L1,L2)))
print('L1 == L3 (before): %s' % str(np.array_equal(L1,L3)))
print('L1 == L4 (before): %s' % str(np.array_equal(L1,L4)))
print('L1 == L5 (before): %s' % str(np.array_equal(L1,L5)))
print('L2 == L3 (before): %s' % str(np.array_equal(L2,L3)))
print('L2 == L4 (before): %s' % str(np.array_equal(L2,L4)))
print('L2 == L5 (before): %s' % str(np.array_equal(L2,L5)))
print('L3 == L4 (before): %s' % str(np.array_equal(L3,L4)))
print('L3 == L5 (before): %s' % str(np.array_equal(L3,L5)))
print('L4 == L5 (before): %s' % str(np.array_equal(L4,L5)))
print('')
# Normally, only np.random.seed(5) should set the random generator for both
# numpy and scipy functions, but just in case, I've included both here
np.random.seed(5)
scipy.random.seed(5)
[eigval,U] = eigh(L,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval1,U1] = eigh(L1,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval2,U2] = eigh(L2,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval3,U3] = eigh(L3,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval4,U4] = eigh(L4,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval5,U5] = eigh(L5,eigvals=(0,1))
# Make sure the inputs still match each other after the function
print('L == L1 (after): %s' % str(np.array_equal(L,L1)))
print('L == L2 (after): %s' % str(np.array_equal(L,L2)))
print('L == L3 (after): %s' % str(np.array_equal(L,L3)))
print('L == L4 (after): %s' % str(np.array_equal(L,L4)))
print('L == L5 (after): %s' % str(np.array_equal(L,L5)))
print('L1 == L2 (after): %s' % str(np.array_equal(L1,L2)))
print('L1 == L3 (after): %s' % str(np.array_equal(L1,L3)))
print('L1 == L4 (after): %s' % str(np.array_equal(L1,L4)))
print('L1 == L5 (after): %s' % str(np.array_equal(L1,L5)))
print('L2 == L3 (after): %s' % str(np.array_equal(L2,L3)))
print('L2 == L4 (after): %s' % str(np.array_equal(L2,L4)))
print('L2 == L5 (after): %s' % str(np.array_equal(L2,L5)))
print('L3 == L4 (after): %s' % str(np.array_equal(L3,L4)))
print('L3 == L5 (after): %s' % str(np.array_equal(L3,L5)))
print('L4 == L5 (after): %s' % str(np.array_equal(L4,L5)))
print('')
# Check to see if the outputs match each other
print('U == U1: %s' % str(np.array_equal(U,U1)))
print('U == U2: %s' % str(np.array_equal(U,U2)))
print('U == U3: %s' % str(np.array_equal(U,U3)))
print('U == U4: %s' % str(np.array_equal(U,U4)))
print('U == U5: %s' % str(np.array_equal(U,U5)))
print('U1 == U2: %s' % str(np.array_equal(U1,U2)))
print('U1 == U3: %s' % str(np.array_equal(U1,U3)))
print('U1 == U4: %s' % str(np.array_equal(U1,U4)))
print('U1 == U5: %s' % str(np.array_equal(U1,U5)))
print('U2 == U3: %s' % str(np.array_equal(U2,U3)))
print('U2 == U4: %s' % str(np.array_equal(U2,U4)))
print('U2 == U5: %s' % str(np.array_equal(U2,U5)))
print('U3 == U4: %s' % str(np.array_equal(U3,U4)))
print('U3 == U5: %s' % str(np.array_equal(U3,U5)))
print('U4 == U5: %s' % str(np.array_equal(U4,U5)))
print('')
print('eigval == eigval1: %s' % str(np.array_equal(eigval,eigval1)))
print('eigval == eigval2: %s' % str(np.array_equal(eigval,eigval2)))
print('eigval == eigval3: %s' % str(np.array_equal(eigval,eigval3)))
print('eigval == eigval4: %s' % str(np.array_equal(eigval,eigval4)))
print('eigval == eigval5: %s' % str(np.array_equal(eigval,eigval5)))
print('eigval1 == eigval2: %s' % str(np.array_equal(eigval1,eigval2)))
print('eigval1 == eigval3: %s' % str(np.array_equal(eigval1,eigval3)))
print('eigval1 == eigval4: %s' % str(np.array_equal(eigval1,eigval4)))
print('eigval1 == eigval5: %s' % str(np.array_equal(eigval1,eigval5)))
print('eigval2 == eigval3: %s' % str(np.array_equal(eigval2,eigval3)))
print('eigval2 == eigval4: %s' % str(np.array_equal(eigval2,eigval4)))
print('eigval2 == eigval5: %s' % str(np.array_equal(eigval2,eigval5)))
print('eigval3 == eigval4: %s' % str(np.array_equal(eigval3,eigval4)))
print('eigval3 == eigval5: %s' % str(np.array_equal(eigval3,eigval5)))
print('eigval4 == eigval5: %s' % str(np.array_equal(eigval4,eigval5)))
print('')
print('#'*50)
这是输出：
##################################################
L symmetric: True
L == L1 (before): True
L == L2 (before): True
L == L3 (before): True
L == L4 (before): True
L == L5 (before): True
L1 == L2 (before): True
L1 == L3 (before): True
L1 == L4 (before): True
L1 == L5 (before): True
L2 == L3 (before): True
L2 == L4 (before): True
L2 == L5 (before): True
L3 == L4 (before): True
L3 == L5 (before): True
L4 == L5 (before): True
L == L1 (after): True
L == L2 (after): True
L == L3 (after): True
L == L4 (after): True
L == L5 (after): True
L1 == L2 (after): True
L1 == L3 (after): True
L1 == L4 (after): True
L1 == L5 (after): True
L2 == L3 (after): True
L2 == L4 (after): True
L2 == L5 (after): True
L3 == L4 (after): True
L3 == L5 (after): True
L4 == L5 (after): True
U == U1: False
U == U2: False
U == U3: False
U == U4: False
U == U5: False
U1 == U2: False
U1 == U3: False
U1 == U4: False
U1 == U5: False
U2 == U3: True
U2 == U4: True
U2 == U5: True
U3 == U4: True
U3 == U5: True
U4 == U5: True
eigval == eigval1: True
eigval == eigval2: True
eigval == eigval3: True
eigval == eigval4: True
eigval == eigval5: True
eigval1 == eigval2: True
eigval1 == eigval3: True
eigval1 == eigval4: True
eigval1 == eigval5: True
eigval2 == eigval3: True
eigval2 == eigval4: True
eigval2 == eigval5: True
eigval3 == eigval4: True
eigval3 == eigval5: True
eigval4 == eigval5: True
##################################################
如何在特征分解上获得一致的、确定性的结果？
因此，根据我收到的评论，我决定尝试更新我的 numpy 和 scipy 库。我将 numpy 从 1.7.1 更新到 1.9.2，将 scipy 从 0.12.0 更新到 0.15.1。这似乎已经成功了！我想这是由于这些库的更高版本中修复了一个错误。
谢谢大家的评论！如果不是因为这个，我可能不会想到更新。
您期望浮点计算总是对相同的输入产生相同的结果实际上是错误的 --- 即使在没有线程的单台计算机上也是如此。如果不采取某些预防措施（会影响性能），高性能计算通常无法实现浮点重现性，您在编写代码时需要考虑到这一点。
这有几个原因，其中之一是优化编译器。
有关详细信息，请参见示例 https://software.intel.com/sites/default/files/article/164389/fp-consistency-102511.pdf。
我从 scipy.linalg.eigh 函数中得到不一致的结果（如果我尝试使用 numpy.linalg.eigh 也是如此）。即使我使用相同的输入，我得到的值也是不同的。为了以防万一，我尝试将随机种子设置为一个常量，但还是不行。
这是我的测试代码：
print('#'*50)
import numpy as np
import scipy
from scipy.linalg import eigh
# Make a random symmetric matrix
m_size = 1500
A = scipy.random.rand(m_size, m_size)
L = np.triu(A) + np.triu(A,1).T
# Verify that matrix is symmetric
print('L symmetric: %s' % np.array_equal(L,L.T))
print('')
# Just to make sure that changes to one L won't effect changes to another
# and that we are starting with the exact same input.
import copy
L1 = copy.deepcopy(L)
L2 = copy.deepcopy(L)
L3 = copy.deepcopy(L)
L4 = copy.deepcopy(L)
L5 = copy.deepcopy(L)
# Sanity check that inputs are the same
print('L == L1 (before): %s' % str(np.array_equal(L,L1)))
print('L == L2 (before): %s' % str(np.array_equal(L,L2)))
print('L == L3 (before): %s' % str(np.array_equal(L,L3)))
print('L == L4 (before): %s' % str(np.array_equal(L,L4)))
print('L == L5 (before): %s' % str(np.array_equal(L,L5)))
print('L1 == L2 (before): %s' % str(np.array_equal(L1,L2)))
print('L1 == L3 (before): %s' % str(np.array_equal(L1,L3)))
print('L1 == L4 (before): %s' % str(np.array_equal(L1,L4)))
print('L1 == L5 (before): %s' % str(np.array_equal(L1,L5)))
print('L2 == L3 (before): %s' % str(np.array_equal(L2,L3)))
print('L2 == L4 (before): %s' % str(np.array_equal(L2,L4)))
print('L2 == L5 (before): %s' % str(np.array_equal(L2,L5)))
print('L3 == L4 (before): %s' % str(np.array_equal(L3,L4)))
print('L3 == L5 (before): %s' % str(np.array_equal(L3,L5)))
print('L4 == L5 (before): %s' % str(np.array_equal(L4,L5)))
print('')
# Normally, only np.random.seed(5) should set the random generator for both
# numpy and scipy functions, but just in case, I've included both here
np.random.seed(5)
scipy.random.seed(5)
[eigval,U] = eigh(L,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval1,U1] = eigh(L1,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval2,U2] = eigh(L2,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval3,U3] = eigh(L3,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval4,U4] = eigh(L4,eigvals=(0,1))
np.random.seed(5)
scipy.random.seed(5)
[eigval5,U5] = eigh(L5,eigvals=(0,1))
# Make sure the inputs still match each other after the function
print('L == L1 (after): %s' % str(np.array_equal(L,L1)))
print('L == L2 (after): %s' % str(np.array_equal(L,L2)))
print('L == L3 (after): %s' % str(np.array_equal(L,L3)))
print('L == L4 (after): %s' % str(np.array_equal(L,L4)))
print('L == L5 (after): %s' % str(np.array_equal(L,L5)))
print('L1 == L2 (after): %s' % str(np.array_equal(L1,L2)))
print('L1 == L3 (after): %s' % str(np.array_equal(L1,L3)))
print('L1 == L4 (after): %s' % str(np.array_equal(L1,L4)))
print('L1 == L5 (after): %s' % str(np.array_equal(L1,L5)))
print('L2 == L3 (after): %s' % str(np.array_equal(L2,L3)))
print('L2 == L4 (after): %s' % str(np.array_equal(L2,L4)))
print('L2 == L5 (after): %s' % str(np.array_equal(L2,L5)))
print('L3 == L4 (after): %s' % str(np.array_equal(L3,L4)))
print('L3 == L5 (after): %s' % str(np.array_equal(L3,L5)))
print('L4 == L5 (after): %s' % str(np.array_equal(L4,L5)))
print('')
# Check to see if the outputs match each other
print('U == U1: %s' % str(np.array_equal(U,U1)))
print('U == U2: %s' % str(np.array_equal(U,U2)))
print('U == U3: %s' % str(np.array_equal(U,U3)))
print('U == U4: %s' % str(np.array_equal(U,U4)))
print('U == U5: %s' % str(np.array_equal(U,U5)))
print('U1 == U2: %s' % str(np.array_equal(U1,U2)))
print('U1 == U3: %s' % str(np.array_equal(U1,U3)))
print('U1 == U4: %s' % str(np.array_equal(U1,U4)))
print('U1 == U5: %s' % str(np.array_equal(U1,U5)))
print('U2 == U3: %s' % str(np.array_equal(U2,U3)))
print('U2 == U4: %s' % str(np.array_equal(U2,U4)))
print('U2 == U5: %s' % str(np.array_equal(U2,U5)))
print('U3 == U4: %s' % str(np.array_equal(U3,U4)))
print('U3 == U5: %s' % str(np.array_equal(U3,U5)))
print('U4 == U5: %s' % str(np.array_equal(U4,U5)))
print('')
print('eigval == eigval1: %s' % str(np.array_equal(eigval,eigval1)))
print('eigval == eigval2: %s' % str(np.array_equal(eigval,eigval2)))
print('eigval == eigval3: %s' % str(np.array_equal(eigval,eigval3)))
print('eigval == eigval4: %s' % str(np.array_equal(eigval,eigval4)))
print('eigval == eigval5: %s' % str(np.array_equal(eigval,eigval5)))
print('eigval1 == eigval2: %s' % str(np.array_equal(eigval1,eigval2)))
print('eigval1 == eigval3: %s' % str(np.array_equal(eigval1,eigval3)))
print('eigval1 == eigval4: %s' % str(np.array_equal(eigval1,eigval4)))
print('eigval1 == eigval5: %s' % str(np.array_equal(eigval1,eigval5)))
print('eigval2 == eigval3: %s' % str(np.array_equal(eigval2,eigval3)))
print('eigval2 == eigval4: %s' % str(np.array_equal(eigval2,eigval4)))
print('eigval2 == eigval5: %s' % str(np.array_equal(eigval2,eigval5)))
print('eigval3 == eigval4: %s' % str(np.array_equal(eigval3,eigval4)))
print('eigval3 == eigval5: %s' % str(np.array_equal(eigval3,eigval5)))
print('eigval4 == eigval5: %s' % str(np.array_equal(eigval4,eigval5)))
print('')
print('#'*50)
这是输出：
##################################################
L symmetric: True
L == L1 (before): True
L == L2 (before): True
L == L3 (before): True
L == L4 (before): True
L == L5 (before): True
L1 == L2 (before): True
L1 == L3 (before): True
L1 == L4 (before): True
L1 == L5 (before): True
L2 == L3 (before): True
L2 == L4 (before): True
L2 == L5 (before): True
L3 == L4 (before): True
L3 == L5 (before): True
L4 == L5 (before): True
L == L1 (after): True
L == L2 (after): True
L == L3 (after): True
L == L4 (after): True
L == L5 (after): True
L1 == L2 (after): True
L1 == L3 (after): True
L1 == L4 (after): True
L1 == L5 (after): True
L2 == L3 (after): True
L2 == L4 (after): True
L2 == L5 (after): True
L3 == L4 (after): True
L3 == L5 (after): True
L4 == L5 (after): True
U == U1: False
U == U2: False
U == U3: False
U == U4: False
U == U5: False
U1 == U2: False
U1 == U3: False
U1 == U4: False
U1 == U5: False
U2 == U3: True
U2 == U4: True
U2 == U5: True
U3 == U4: True
U3 == U5: True
U4 == U5: True
eigval == eigval1: True
eigval == eigval2: True
eigval == eigval3: True
eigval == eigval4: True
eigval == eigval5: True
eigval1 == eigval2: True
eigval1 == eigval3: True
eigval1 == eigval4: True
eigval1 == eigval5: True
eigval2 == eigval3: True
eigval2 == eigval4: True
eigval2 == eigval5: True
eigval3 == eigval4: True
eigval3 == eigval5: True
eigval4 == eigval5: True
##################################################
如何在特征分解上获得一致的、确定性的结果？
因此，根据我收到的评论，我决定尝试更新我的 numpy 和 scipy 库。我将 numpy 从 1.7.1 更新到 1.9.2，将 scipy 从 0.12.0 更新到 0.15.1。这似乎已经成功了！我想这是由于这些库的更高版本中修复了一个错误。
谢谢大家的评论！如果不是因为这个，我可能不会想到更新。
您期望浮点计算总是对相同的输入产生相同的结果实际上是错误的 --- 即使在没有线程的单台计算机上也是如此。如果不采取某些预防措施（会影响性能），高性能计算通常无法实现浮点重现性，您在编写代码时需要考虑到这一点。
这有几个原因，其中之一是优化编译器。 有关详细信息，请参见示例 https://software.intel.com/sites/default/files/article/164389/fp-consistency-102511.pdf。