Детерминированный python script ведет себя недетерминированным способом
У меня есть script, который не использует рандомизацию, которая дает мне разные ответы при запуске. Я ожидаю, что ответ будет таким же, каждый раз, когда я запускаю script. Проблема, по-видимому, возникает только для определенных (плохо обусловленных) входных данных.
Отрывок из алгоритма используется для вычисления определенного типа контроллера для линейной системы и в основном состоит из выполнения линейной алгебры (матричные инверсии, уравнение Риккати, собственные значения).
Очевидно, это серьезное беспокойство для меня, так как теперь я не могу доверять своему коду, чтобы дать мне правильные результаты. Я знаю, что результат может быть неправильным для плохо подготовленных данных, но я ожидаю, что это будет неправильно. Почему ответ на моей машине Windows не всегда одинаковый? Почему машины Linux и Windows не дают одинаковых результатов?
Я использую Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win 32, с Numpy версии 1.8.2 и Scipy 0.14.0. (Windows 8, 64 бит).
Код ниже. Я также попытался запустить код на двух машинах Linux, и там script всегда дает тот же ответ (но машины дали разные ответы). Один из них запускал Python 2.7.8, с Numpy 1.8.2 и Scipy 0.14.0. Второй выполнял Python 2.7.3 с Numpy 1.6.1 и Scipy 0.12.0.
Я решаю уравнение Риккати три раза, а затем печатаю ответы. Я ожидаю такой же ответ каждый раз, вместо этого получаю последовательность "1.75305103767e-09; 3.25501787302e-07; 3.25501787302e-07'.
import numpy as np
import scipy.linalg
matrix = np.matrix
A = matrix([[ 0.00000000e+00, 2.96156260e+01, 0.00000000e+00,
-1.00000000e+00],
[ -2.96156260e+01, -6.77626358e-21, 1.00000000e+00,
-2.11758237e-22],
[ 0.00000000e+00, 0.00000000e+00, 2.06196064e+00,
5.59422224e+01],
[ 0.00000000e+00, 0.00000000e+00, 2.12407340e+01,
-2.06195974e+00]])
B = matrix([[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ -342.35401351, -14204.86532216, 31.22469724],
[ 1390.44997337, 342.33745324, -126.81720597]])
Q = matrix([[ 5.00000001, 0. , 0. , 0. ],
[ 0. , 5.00000001, 0. , 0. ],
[ 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. ]])
R = matrix([[ -3.75632852e+04, -0.00000000e+00, 0.00000000e+00],
[ -0.00000000e+00, -3.75632852e+04, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 4.00000000e+00]])
counter = 0
while counter < 3:
counter +=1
X = scipy.linalg.solve_continuous_are(A, B, Q, R)
print(-3449.15531628 - X[0,0])
Моя конфигурация numpy ниже print np.show_config()
lapack_opt_info:
libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md', 'mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md']
library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32']
define_macros = [('SCIPY_MKL_H', None)]
include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include']
blas_opt_info:
libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md']
library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32']
define_macros = [('SCIPY_MKL_H', None)]
include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include']
openblas_info:
NOT AVAILABLE
lapack_mkl_info:
libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md', 'mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md']
library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32']
define_macros = [('SCIPY_MKL_H', None)]
include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include']
blas_mkl_info:
libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md']
library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32']
define_macros = [('SCIPY_MKL_H', None)]
include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include']
mkl_info:
libraries = ['mkl_blas95', 'mkl_lapack95', 'mkl_intel_c', 'mkl_intel_thread', 'mkl_core', 'libiomp5md']
library_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/lib/ia32', 'C:/Program Files (x86)/Intel/Composer XE 2013 SP1/compiler/lib/ia32']
define_macros = [('SCIPY_MKL_H', None)]
include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include']
None
(редактирует, чтобы обрезать вопрос)
Ответы
Ответ 1
В общем, библиотеки linalg в Windows дают разные ответы на разные прогоны на уровне точности машины. Я никогда не слышал объяснений, почему это происходит только или в основном на Windows.
Если ваша матрица плохо кондиционирована, то inv будет в значительной степени численным шумом. В Windows шум не всегда одинаковый в последовательных прогонах, в других операционных системах шум может быть всегда одним и тем же, но может различаться в зависимости от деталей библиотеки линейной алгебры, от параметров потоков, использования кеша и т.д.
Я видел и отправлял в scipy список рассылки несколько примеров для этого в Windows, я использовал в основном официальные 32-битные двоичные файлы с ATLAS BLAS/LAPACK.
Единственное решение состоит в том, чтобы сделать результат вашего расчета не так сильно зависящим от точности точности с плавающей запятой и числовым шумом, например, упорядочить матрицу обратным, использовать обобщенный инверсный, pinv, reparameterize или аналогичный.
Ответ 2
Как pv, отмеченный в комментариях к answer33700 answer, предыдущая формулировка решателей Riccati были, хотя и теоретически правильными, не численно устойчивыми. Эта проблема исправлена в версии разработки SciPy, а обобщающие уравнения Riccati также поддерживают обобщения.
Репликаторы Riccati обновляются, и результирующие решатели будут доступны с версии 0.19 и далее. Здесь вы можете проверить .
Если, используя данный пример в вопросе, я заменю последний цикл на
for _ in range(5):
x = scipy.linalg.solve_continuous_are(A, B, Q, R)
Res = [email protected] + [email protected] + q - [email protected]@ np.linalg.solve(r,b.T)@ x
print(Res)
Я получаю (окна 10, py3.5.2)
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06]
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05]
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06]
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]]
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06]
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05]
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06]
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]]
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06]
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05]
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06]
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]]
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06]
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05]
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06]
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]]
[[ 2.32314924e-05 -2.55086270e-05 -7.66709854e-06 -9.01878229e-06]
[ -2.62447211e-05 2.61182140e-05 8.27328768e-06 1.00345896e-05]
[ -7.92257197e-06 8.57094892e-06 2.50908488e-06 3.05714639e-06]
[ -9.51046241e-06 9.80847472e-06 3.13103374e-06 3.60747799e-06]]
Для справки, возвращаемое решение
array([[-3449.15531305, 4097.1707738 , 1142.52971904, 1566.51333847],
[ 4097.1707738 , -4863.70533241, -1356.66524959, -1860.15980716],
[ 1142.52971904, -1356.66524959, -378.32586814, -518.71965102],
[ 1566.51333847, -1860.15980716, -518.71965102, -711.21062988]])
