Непоследовательные результаты для функции Numpy/Scipy eigh

Question

Я получаю несогласованные результаты из функции 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 

################################################## 

Как я могу получить последовательные, детерминированных содержат результаты на eigendecomposition?

Answer 1

Ваше ожидание того, что вычисления с плавающей запятой всегда дают одинаковый результат для одних и тех же входов, на самом деле ошибочны - даже на одном компьютере без потоков. Воспроизводимость с плавающей точкой вообще не достигается при высокопроизводительных вычислениях без определенных мер предосторожности (которые имеют влияние на производительность), которые необходимо учитывать при написании кода.

Существует несколько причин для этого, один из которых является оптимизацией компиляторов. См., Например, https://software.intel.com/sites/default/files/article/164389/fp-consistency-102511.pdf.

Answer 2

Итак, основываясь на полученных комментариях, я решил попробовать обновить свои библиотеки numpy и scipy. Я обновил numpy от 1.7.1 до 1.9.2 и scipy от 0.12.0 до 0.15.1. И это похоже на трюк! Полагаю, это связано с ошибкой, которая исправлена ​​в более поздних версиях этих библиотек.

Спасибо за комментарии, ребята! Возможно, я бы не подумал обновить, если бы не это.