I am getting inconsistent results from the scipy.linalg.eigh function (and also if I try using numpy.linalg.eigh). Even when I use the same input, the values that I get out are different. Just in case, I tried setting the random seed to a constant, but still no luck.
Here is my test code:
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)
And here is the output:
##################################################
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
##################################################
How can I get consistent, deterministic results on an eigendecomposition?
Your expectation that floating point calculations always yield the same result for same inputs is actually wrong --- even on a single computer without threads. Floating point reproducibility in general not achieved in high-performance computing without certain precautions (which have performance impact), which you need to take into account when writing code.
There are several reasons for this, one of them being optimizing compilers. See for example https://software.intel.com/sites/default/files/article/164389/fp-consistency-102511.pdf for details.