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)


Answer: 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.

