确定性 python 脚本以非确定性方式运行答案

【问题标题】：Deterministic python script behaves in non-deterministic way确定性 python 脚本以非确定性方式运行
【发布时间】：2015-02-19 10:28:06
【问题描述】：

我有一个不使用随机化的脚本，当我运行它时会给我不同的答案。我希望每次运行脚本时答案都是一样的。这个问题似乎只发生在某些（病态的）输入数据上。

sn-p 来自为线性系统计算特定类型控制器的算法，它主要由线性代数（矩阵求逆、Riccati 方程、特征值）组成。

显然，这对我来说是一个主要的担忧，因为我现在不能相信我的代码会给我正确的结果。我知道对于条件不佳的数据，结果可能是错误的，但我希望始终是错误的。为什么我的 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。第二个是使用 Numpy 1.6.1 和 Scipy 0.12.0 运行 Python 2.7.3。

我求解 Riccati 方程 3 次，然后打印答案。我每次都希望得到相同的答案，而是得到序列 '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：库= ['mkl_blas95'，'mkl_lapack95'，'mkl_intel_c'，'mkl_intel_thread'，'mkl_core'，'libiomp5md'，'mkl_blas95'，'mkl_lapack95'，'mkl_intel_c'，'mkl_intel_thread'，'d'libiomp5'， ] 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' ] 定义宏 = [（'SCIPY_MKL_H'，无）] include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] blas_opt_info：库 = ['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' ] 定义宏 = [（'SCIPY_MKL_H'，无）] include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] openblas_info：无法使用 lapack_mkl_info：库= ['mkl_blas95'，'mkl_lapack95'，'mkl_intel_c'，'mkl_intel_thread'，'mkl_core'，'libiomp5md'，'mkl_blas95'，'mkl_lapack95'，'mkl_intel_c'，'mkl_intel_thread'，'d'libiomp5'， ] 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' ] 定义宏 = [（'SCIPY_MKL_H'，无）] include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] blas_mkl_info：库 = ['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' ] 定义宏 = [（'SCIPY_MKL_H'，无）] include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] mkl_info：库 = ['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', 无)] include_dirs = ['c:/Program Files (x86)/Intel/Composer XE 2013 SP1/mkl/include'] 没有

（编辑以减少问题）

【问题讨论】：

如果您将日志语句添加到代码中以报告各个阶段的关键变量的值，这将有很大帮助。然后你可以在不同机器生成的日志文件上运行diff，看看代码在哪里出现了分歧。
我正要建议@unutbu 的提议。这应该有助于缩小问题范围。
我们的目标是将代码筛选成最简单的具体示例，在不同的机器上产生不同的值。
您的脚本在我的系统上始终给出 -305367.072356（Windows 7、Python 2.7.8 (Anaconda 2.0.1)、Numpy：1.8.1、Scipy：0.14.0）。很奇怪。
一个疯狂的猜测 - 底层实现是并行/多线程的吗？如果是这样，可能是一个并发错误 - 可以解释不确定的行为（但可能会看到与不同机器上的内核数量相关）？

标签： python windows numpy scipy non-deterministic

【解决方案1】：

一般来说，Windows 上的 linalg 库在机器精度级别的不同运行中给出不同的答案。我从来没有听说过为什么这种情况只发生或主要发生在 Windows 上。

如果您的矩阵是病态的，那么 inv 将主要是数值噪声。在 Windows 上，连续运行的噪声并不总是相同，在其他操作系统上，噪声可能总是相同，但可能会有所不同，具体取决于线性代数库的细节、线程选项、缓存使用情况等。

我在 scipy 邮件列表上看到并发布了几个 Windows 上的示例，我主要使用带有 ATLAS BLAS/LAPACK 的官方 32 位二进制文件。

唯一的解决方案是使您的计算结果不太依赖于浮点精度问题和数值噪声，例如正则化矩阵逆，使用广义逆，pinv，重新参数化或类似的。

【讨论】：

部分数值噪音来自内存中的数据对齐，即是否恰好位于适合cpu的sse等指令的位置---这两种情况的操作顺序不同，所以浮点舍入误差是不同的。 win32 标准内存分配器 afaik 倾向于返回平均比 linux 分配器对齐更少的内存块，因此那里有更多的抖动。如果您的结果很大程度上取决于 FP 舍入误差的走向，那么您的问题表述在数值上是不稳定的。

【解决方案2】：

正如pv 在comments to user333700's answer 中指出的那样，Riccati 求解器的先前公式虽然在理论上是正确的，但在数值上并不稳定。此问题已在 SciPy 的开发版本中得到修复，并且求解器也支持广义 Riccati 方程。

Riccati 求解器已更新，生成的求解器将从 0.19 及更高版本开始提供。您可以查看development branch docs here。

如果，使用问题中的给定示例，我将最后一个循环替换为

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)

我得到（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]]

供参考，返回的解决方案是

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

【讨论】：