Deterministic python script behaves in non-deterministic way

Question

I have a script which uses no randomisation that gives me different answers when I run it. I expect the answer to be the same, every time I run the script. The problem appears to only happen for certain (ill-conditioned) input data.

The snippet comes from an algorithm to compute a specific type of controller for a linear system, and it mostly consists of doing linear algebra (matrix inversions, Riccati equation, eigenvalues).

Obviously, this is a major worry for me, as I now cannot trust my code to give me the right results. I know the result can be wrong for poorly conditioned data, but I expect consistently wrong. Why is the answer not always the same on my Windows machine? Why do the Linux & Windows machine not give the same results?

I'm using Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win 32, with Numpy version 1.8.2 and Scipy 0.14.0. (Windows 8, 64bit).

The code is below. I've also tried running the code on two Linux machines, and there the script always gives the same answer (but the machines gave differing answers). One was running Python 2.7.8, with Numpy 1.8.2 and Scipy 0.14.0. The second was running Python 2.7.3 with Numpy 1.6.1 and Scipy 0.12.0.

I solve the Riccati equation three times, and then print the answers. I expect the same answer every time, instead I get the sequence '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])

My numpy config is as below 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

(edits to trim the question down)

Did you seed the generator? http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html — NPE, Dec 20 '14 at 20:01
I have added `np.random.seed(0)` as the first line of code (just after the inputs), and it makes no difference. Will update question. — mwmwm, Dec 20 '14 at 20:06
Apologies, I somehow misread your question as stating that the code *does* use randomisation. Oops. — NPE, Dec 20 '14 at 20:06
It would help greatly if you add logging statements to your code which report the value of key variables at various stages. Then you could run `diff` on the log files generated by the different machines and see where in the code things diverge. — unutbu, Dec 20 '14 at 20:10
I was just about to suggest what @unutbu is proposing. That should help narrow the problem down. — NPE, Dec 20 '14 at 20:12
The goal is to winnow the code down to the simplest concrete example which produces different values on different machines. — unutbu, Dec 20 '14 at 20:16
Your script consistently gives -305367.072356 on my system (Windows 7, Python 2.7.8 (Anaconda 2.0.1), Numpy: 1.8.1, Scipy: 0.14.0). Weird. — Roberto, Dec 20 '14 at 21:10
A wild guess - are the underlying implementations parallel/multi-threaded? If so, could be a concurrency bug - would explain non-deterministic behaviour (but might see a correletion with the number of cores on different machines)? — DNA, Dec 20 '14 at 21:31
I also can't reproduce the non-deterministic behavior. Using numpy 1.8.2 / scipy 0.13.3 / Python 2.7.6 running on Ubuntu 14.04 I consistently get -255097.48248. Running newer dev versions (numpy 1.10.0.dev-8bcb756 and scipy 0.15.0.dev-aacfbfb) I consistently get -3608625.51589. Can you tell us which BLAS library you're using on your Windows machine? What's the output of `np.show_config()`? — ali_m, Dec 20 '14 at 23:29
I've added a smaller example, focused on the Riccati equation where the mischief starts. I've also added my numpy config output. — mwmwm, Dec 21 '14 at 06:23
@mwmwm I would suggest trimming out the old, larger version. — Veedrac, Dec 21 '14 at 19:49

score 8 · Accepted Answer · answered Dec 21 '14 at 05:17

In general, linalg libraries on Windows give different answers on different runs at machine precision level. I never heard of an explanation why this happens only or mainly on Windows.

If your matrix is ill conditioned, then the inv will be largely numerical noise. On Windows the noise is not always the same in consecutive runs, on other operating systems the noise might be always the same but can differ depending on the details of the linear algebra library, on threading options, cache usage and so on.

I've seen on and posted to the scipy mailing list several examples for this on Windows, I was using mostly the official 32 bit binaries with ATLAS BLAS/LAPACK.

The only solution is to make the outcome of your calculation not depend so much on floating point precision issues and numerical noise, for example regularize the matrix inverse, use generalized inverse, pinv, reparameterize or similar.

Part of the numerical noise comes from data alignment in memory, ie whether it happens to lie directly in a position suitable for cpu's sse etc. instructions or not --- the two cases do operations in different order, so floating point rounding error is different. win32 standard memory allocator afaik tends to return blocks of memory that are on average less aligned than linux allocators, so there is more jitter there. If your results depend significantly on which way FP rounding error goes, your problem formulation is numerically unstable. — pv., Dec 29 '14 at 13:46

score 2 · Answer 2 · edited May 23 '17 at 11:53

As pv noted in the comments to user333700's answer, the previous formulation of the Riccati solvers were, though being theoretically correct, not numerically stable. This issue is fixed on the development version of SciPy and the solvers support generalized Riccati equations too.

The Riccati solvers are updated and resulting solvers will be available from version 0.19 and onwards. You can check the development branch docs here.

If, using the given example in the question I replace the last loop with

for _ in range(5):
    x = scipy.linalg.solve_continuous_are(A, B, Q, R)
    Res = x@a + a.T@x + q - x@b@ np.linalg.solve(r,b.T)@ x
    print(Res)

I get (windows 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]]

For reference, the solution returned is

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]])

Deterministic python script behaves in non-deterministic way

2 Answers2

Linked