python - para - instalar pip en windows
El script de python determinista se comporta de manera no determinista (2)
Tengo un script que no usa asignación al azar que me da respuestas diferentes cuando lo ejecuto. Espero que la respuesta sea la misma, cada vez que ejecuto el script. El problema parece ocurrir solo para ciertos datos de entrada (mal acondicionados).
El fragmento proviene de un algoritmo para calcular un tipo específico de controlador para un sistema lineal, y consiste principalmente en hacer álgebra lineal (inversiones de matriz, ecuación de Riccati, valores propios).
Obviamente, esta es una gran preocupación para mí, ya que ahora no puedo confiar en mi código para darme los resultados correctos. Sé que el resultado puede ser incorrecto para los datos mal condicionados, pero espero que siempre esté mal. ¿Por qué la respuesta no siempre es la misma en mi máquina con Windows? ¿Por qué la máquina de Linux y Windows no da los mismos resultados?
Estoy usando Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win 32
, con Numpy versión 1.8.2 y Scipy 0.14.0. (Windows 8, 64 bits).
El código está abajo. También intenté ejecutar el código en dos máquinas Linux, y allí el script siempre da la misma respuesta (pero las máquinas dieron respuestas diferentes). Uno ejecutaba Python 2.7.8, con Numpy 1.8.2 y Scipy 0.14.0. El segundo estaba ejecutando Python 2.7.3 con Numpy 1.6.1 y Scipy 0.12.0.
Resuelvo la ecuación de Riccati tres veces y luego imprimo las respuestas. Espero la misma respuesta todas las veces, en vez de eso obtengo la secuencia ''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])
Mi configuración numpy es la siguiente, 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
(ediciones para recortar la pregunta)
Como pv notó en los comentarios a la respuesta del usuario333700 , la formulación anterior de los solucionadores de Riccati eran, aunque teóricamente correctos, no numéricamente estables. Este problema está fijado en la versión de desarrollo de SciPy y los solucionadores también soportan ecuaciones de Riccati generalizadas.
Los solucionadores de Riccati se actualizan y los solucionadores resultantes estarán disponibles desde la versión 0.19 en adelante. Puede consultar los documentos de la sucursal de desarrollo aquí .
Si, usando el ejemplo dado en la pregunta, reemplazo el último ciclo con
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)
Me sale (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]]
Como referencia, la solución devuelta es
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]])
En general, las bibliotecas de linalg en Windows dan diferentes respuestas en diferentes ejecuciones al nivel de precisión de la máquina. Nunca escuché una explicación de por qué esto sucede solo o principalmente en Windows.
Si su matriz está mal condicionada, entonces el inv será en gran parte ruido numérico. En Windows, el ruido no siempre es el mismo en ejecuciones consecutivas, en otros sistemas operativos el ruido puede ser siempre el mismo, pero puede variar según los detalles de la biblioteca de álgebra lineal, las opciones de subprocesamiento, el uso de la memoria caché, etc.
He visto y publicado en la lista de correo scipy varios ejemplos para esto en Windows, estaba usando principalmente los binarios oficiales de 32 bits con ATLAS BLAS / LAPACK.
La única solución es hacer que el resultado de su cálculo no dependa tanto de cuestiones de precisión de coma flotante y ruido numérico, por ejemplo, regularizar la matriz inversa, usar inversión inversa generalizada, pinv, reparametrizar o similar.