2014-12-20 2 views
12

У меня есть сценарий, который не использует рандомизацию, которая дает мне разные ответы при запуске. Я ожидаю, что ответ будет таким же, каждый раз, когда я запускаю скрипт. Проблема, по-видимому, возникает только для определенных (плохо обусловленных) входных данных.Детерминированный скрипт python ведет себя недетерминированным способом

Отрывок из алгоритма для вычисления определенного типа контроллера для линейной системы и в основном состоит из выполнения линейной алгебры (матричные инверсии, уравнение Риккати, собственные значения).

Очевидно, это серьезное беспокойство для меня, так как теперь я не могу доверять своему коду, чтобы дать мне правильные результаты. Я знаю, что результат может быть неправильным для плохо подготовленных данных, но я ожидаю, что это будет неправильно. Почему ответ на моей машине 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, и там сценарий всегда дает тот же ответ (но машины дали разные ответы). Один из них запускал 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 

(редактирует подрезать вопрос вниз)

+0

Вы заселяли генератор? http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html – NPE

+0

Я добавил 'np.random.seed (0)' в качестве первой строки кода (сразу после ввода), и это не имеет значения. Будет обновлен вопрос. – mwmwm

+0

Извинения, я как-то неправильно прочитал ваш вопрос, заявив, что код * делает * использует рандомизацию. К сожалению. – NPE

ответ

7

В общем, linalg библиотеки на Windows, дают разные ответы на разных прогонов на уровне точности станка , Я никогда не слышал объяснений, почему это происходит только или в основном на Windows.

Если ваша матрица плохо кондиционирована, то inv будет в значительной степени численным шумом. В Windows шум не всегда одинаковый в последовательных прогонах, в других операционных системах шум может быть всегда одним и тем же, но может различаться в зависимости от деталей библиотеки линейной алгебры, параметров потоков, использования кеша и т. Д.

Я видел и отправлял в scipy список рассылки несколько примеров для этого в Windows, я использовал в основном официальные 32-битные двоичные файлы с ATLAS BLAS/LAPACK.

Единственное решение состоит в том, чтобы сделать результат вашего расчета не так сильно зависящим от точности точности с плавающей запятой и числовым шумом, например, упорядочить матрицу обратным, использовать обобщенный инверсный, pinv, reparameterize или аналогичный.

+3

Часть числового шума возникает из-за выравнивания данных в памяти, т. Е. Является ли оно непосредственно лежащим в позиции, подходящей для инструкций sse и т. Д., Или нет. Эти два случая выполняют операции в другом порядке, поэтому ошибка округления с плавающей запятой другой. В стандартном распределении памяти win32 afaik имеет тенденцию возвращать блоки памяти, которые в среднем меньше выровнены, чем линейные распределители, поэтому там больше джиттера. Если ваши результаты существенно зависят от того, по каким причинам происходит ошибка округления FP, ваша формулировка проблемы численно неустойчива. –

2

Как указано в comments to user333700's answer, предыдущая формулировка решателей Риккати была, хотя и теоретически правильной, не была численно устойчивой. Эта проблема исправлена ​​в версии разработки SciPy, а обобщающие уравнения Riccati также поддерживают обобщения.

Решители Riccati обновляются, и результирующие решатели будут доступны с версии 0.19 и далее. Вы можете проверить development branch docs here.

Если, используя данный пример в этом вопросе я заменить последнюю петлю с

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]]) 
Смежные вопросы