Deterministic python script behaves in non-deterministic way

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.

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