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 1

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.

Answer 2

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