我有一个脚本,不使用随机化,当我运行它给了我不同的答案。 我希望答案是一样的,每次运行脚本。 这个问题似乎只发生在某些(病态的)input数据上。
片段来自一个algorithm来计算一个线性系统的特定types的控制器,它主要包括做线性代数(matrix求逆,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方程,然后打印答案。 我期望每次都有相同的答案,而不是我得到序列'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])
我的numpyconfiguration如下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','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: 库= ['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: 无法使用 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','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: 库= ['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: 库= ['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'] 没有
(编辑修剪问题)
一般来说,Windows上的linalg库在机器精度级别的不同运行上给出了不同的答案。 我从来没有听说过为什么这只会发生,或主要是在Windows上的解释。
如果你的矩阵有条件的话,那么inv就是数字噪声。 在Windows上,连续运行时噪声并不总是相同,在其他操作系统上,噪声可能总是相同,但可能会因线性代数库的细节,线程选项,缓存使用情况等而有所不同。
我已经看到,并张贴到scipy邮件列表几个这样的例子在Windows上,我主要使用正式的32位二进制文件与ATLAS BLAS / LAPACK。
唯一的解决办法是让计算结果不要太依赖于浮点精度问题和数值噪声,例如调整矩阵的逆矩阵,使用广义逆,pinv,重新参数化或类似的。
正如pv在user333700的回答中所指出的那样,Riccati解算器的先前的表述尽管在理论上是正确的,但在数字上并不稳定。 这个问题固定在SciPy的开发版本上,解算器也支持广义Riccati方程。
Riccati求解器被更新,解析器将从版本0.19开始提供。 你可以在这里查看开发分支文档 。
如果使用问题中给出的例子,我用最后一个循环替换
for _ in range(5): x = scipy.linalg.solve_continuous_are(A, B, Q, R) Res = x@a + aT@x + q - x@b@ np.linalg.solve(r,bT)@ 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]])