Deterministisches Python-Skript verhält sich nicht-deterministisch

Ich habe ein Skript, das keine Randomisierung verwendet, die mir unterschiedliche Antworten gibt, wenn ich es ausführe. Ich erwarte, dass die Antwort immer gleich ist, wenn ich das Skript ausführe. Das Problem scheint nur für bestimmte (schlecht konditionierte) Eingabedaten zu bestehen.Deterministisches Python-Skript verhält sich nicht-deterministisch

Das Snippet kommt von einem Algorithmus, um einen bestimmten Typ von Controller für ein lineares System zu berechnen, und es besteht hauptsächlich aus linearer Algebra (Matrixinversionen, Riccati-Gleichung, Eigenwerte).

Offensichtlich ist dies eine große Sorge für mich, da ich jetzt meinem Code nicht vertrauen kann, um mir die richtigen Ergebnisse zu geben. Ich weiß, dass das Ergebnis für schlecht konditionierte Daten falsch sein kann, aber ich erwarte konsequent falsch. Warum ist die Antwort auf meinem Windows-Rechner nicht immer gleich? Warum gibt der Linux-Rechner & nicht die gleichen Ergebnisse?

Ich benutze Python 2.7.9 (default, Dec 10 2014, 12:24:55) [MSC v.1500 32 bit (Intel)] on win 32, mit Numpy Version 1.8.2 und Scipy 0.14.0. (Windows 8, 64 Bit).

Der Code ist unten. Ich habe auch versucht, den Code auf zwei Linux-Rechnern auszuführen, und dort gibt das Skript immer die gleiche Antwort (aber die Maschinen gaben unterschiedliche Antworten). Eines lief mit Python 2.7.8, mit Numpy 1.8.2 und Scipy 0.14.0. Die zweite war Python 2.7.3 mit Numpy 1.6.1 und Scipy 0.12.0.

Ich löse die Riccati Gleichung dreimal, und drucke dann die Antworten. Ich erwarte immer die gleiche Antwort, stattdessen bekomme ich die Sequenz '1.75305103767e-09; 3.25501787302d-07; 3.25501787302-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])

Meine numpy Config ist als unten 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

(bearbeitet die Frage zu trimmen)

Quelle

2014-12-20 mwmwm

Haben Sie den Generator gesät? http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html – NPE

Ich habe 'np.random.seed (0)' als die erste Zeile des Codes (direkt nach den Eingaben) hinzugefügt), und es macht keinen Unterschied. Wird die Frage aktualisieren. – mwmwm

Entschuldigung, ich habe Ihre Frage irgendwie falsch gelesen und festgestellt, dass der Code * eine Randomisierung verwendet. Hoppla. – NPE

Im Allgemeinen linalg Bibliotheken unter Windows geben unterschiedliche Antworten auf verschiedenen Läufen auf Maschinengenauigkeit Ebene . Ich habe noch nie von einer Erklärung gehört, warum das nur oder hauptsächlich unter Windows passiert.

Wenn Ihre Matrix schlecht konditioniert ist, dann wird das Inv-Signal weitgehend numerisch sein. Unter Windows ist das Rauschen in aufeinanderfolgenden Durchläufen nicht immer gleich, auf anderen Betriebssystemen kann das Rauschen immer gleich sein, kann aber abhängig von den Details der Bibliothek der linearen Algebra, den Threading-Optionen, der Cache-Nutzung usw. abweichen.

Ich habe auf der scipy Mailingliste einige Beispiele dafür auf Windows gesehen und gepostet, ich benutzte hauptsächlich die offiziellen 32 Bit Binärdateien mit ATLAS BLAS/LAPACK. Die einzige Lösung besteht darin, das Ergebnis Ihrer Berechnung nicht so sehr von Gleitkomma-Genauigkeitsproblemen und numerischem Rauschen abhängig zu machen, z. B. die Matrix-Inversität zu regularisieren, verallgemeinerte Inverse, pinv, umparametrisieren oder ähnliches zu verwenden.

Quelle

2014-12-21 05:17:54 user333700

Ein Teil des numerischen Rauschens kommt von der Datenausrichtung im Speicher, dh ob es zufällig direkt an einer Position liegt, die für CPU-Anweisungen usw. geeignet ist oder nicht --- die beiden Fälle führen Operationen in unterschiedlicher Reihenfolge durch, so dass Gleitkomma-Rundungsfehler auftreten anders. win32 standard memory allocator afaik tendiert dazu, Speicherblöcke zurückzugeben, die im Durchschnitt weniger ausgerichtet sind als Linux-Allokatoren, daher gibt es dort mehr Jitter. Wenn Ihre Ergebnisse wesentlich davon abhängen, wie der FP-Rundungsfehler verläuft, ist Ihre Problemformulierung numerisch instabil. –

Wie pv in der comments to user333700's answer notiert, war die vorherige Formulierung der Riccati-Löser, obwohl theoretisch korrekt, nicht numerisch stabil. Dieses Problem wurde in der Entwicklungsversion von SciPy behoben, und die Solver unterstützen auch verallgemeinerte Riccati-Gleichungen.

Die Riccati-Solver werden aktualisiert und die resultierenden Solver sind ab Version 0.19 verfügbar. Sie können die development branch docs here überprüfen.

Wenn das gegebene Beispiel in der Frage ich mit der letzten Schleife ersetzen mit

for _ in range(5): 
    x = scipy.linalg.solve_continuous_are(A, B, Q, R) 
    Res = [email protected] + [email protected] + q - [email protected]@ np.linalg.solve(r,b.T)@ x 
    print(Res)

ich (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]]

Als Referenz zurückgegeben die Lösung

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

Quelle

2016-10-15 15:53:09 percusse

Deterministisches Python-Skript verhält sich nicht-deterministisch

Antwort

Verwandte Themen