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

Question

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

该代码段来自一种计算线性系统特定类型 Controller 的算法，它主要包括线性代数(矩阵求逆、Riccati 方程、特征值)。

显然，这对我来说是一个主要的担忧，因为我现在不能相信我的代码会给我正确的结果。我知道对于条件不佳的数据，结果可能是错误的，但我预计始终是错误的。为什么我的 Windows 机器上的答案并不总是相同？为什么 Linux 和 Windows 机器给出的结果不同？

我在 win 32 上使用 Python 2.7.9(默认，2014 年 12 月 10 日，12:24:55)[MSC v.1500 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 方程，然后打印答案。我每次都希望得到相同的答案，但我得到的是序列 '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:    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 AVAILABLElapack_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

(编辑以减少问题)

Answer 1

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

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

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

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