Matlab/Octave/Numpy numerischen Unterschied

Question

Ich portiere einen Algorithmus, der in Matlab funktioniert, um zu numpy und beobachtete ein seltsames Verhalten. Der relevante Codesegment ist

P = eye(4)*1e20; 
A = [1 -0.015 -0.025 -0.035; 0.015 1 0.035 -0.025; 0.025 -0.035 1 0.015; 0.035 0.025 -0.015 1]; 
V1 = A*(P*A') 
V2 = (A*P)*A' 

Dieser Code, wenn ich mit Matlab laufen, bietet die folgenden Matrizen:

V1 = 1.0021e+20   0 -8.0000e+00   0 
       0 1.0021e+20   0   0 
    -8.0000e+00   0 1.0021e+20   0 
       0   0   0 1.0021e+20 


V2 = 1.0021e+20   0 -8.0000e+00   0 
       0 1.0021e+20   0   0 
    -8.0000e+00   0 1.0021e+20   0 
       0   0   0 1.0021e+20 

Man beachte, daß V1 und V2 gleich sind, wie erwartet.

Wenn derselbe Code in Octave läuft, bietet es:

V1 = 1.0021e+20 4.6172e+01 -1.3800e+02 1.8250e+02 
    -4.6172e+01 1.0021e+20 -1.8258e+02 -1.3800e+02 
    1.3801e+02 1.8239e+02 1.0021e+20 -4.6125e+01 
    -1.8250e+02 1.3800e+02 4.6125e+01 1.0021e+20 

V2 = 1.0021e+20 -4.6172e+01 1.3801e+02 -1.8250e+02 
    4.6172e+01 1.0021e+20 1.8239e+02 1.3800e+02 
    -1.3800e+02 -1.8258e+02 1.0021e+20 4.6125e+01 
    1.8250e+02 -1.3800e+02 -4.6125e+01 1.0021e+20 

In numpy das Segment

from numpy import array, dot, eye 
A = numpy.array([[1, -0.015, -0.025, -0.035],[0.015, 1, 0.035, -0.025],[0.025, -0.035, 1, 0.015],[0.035, 0.025, -0.015, 1]]) 
P = numpy.eye(4)*1e20 
print numpy.dot(A,numpy.dot(P,A.transpose())) 
print numpy.dot(numpy.dot(A,P),A.transpose()) 

die

[[ 1.00207500e+20 4.61718750e+01 -1.37996094e+02 1.82500000e+02] 
[ -4.61718750e+01 1.00207500e+20 -1.82582031e+02 -1.38000000e+02] 
[ 1.38011719e+02 1.82386719e+02 1.00207500e+20 -4.61250000e+01] 
[ -1.82500000e+02 1.38000000e+02 4.61250000e+01 1.00207500e+20]] 
[[ 1.00207500e+20 -4.61718750e+01 1.38011719e+02 -1.82500000e+02] 
[ 4.61718750e+01 1.00207500e+20 1.82386719e+02 1.38000000e+02] 
[ -1.37996094e+02 -1.82582031e+02 1.00207500e+20 4.61250000e+01] 
[ 1.82500000e+02 -1.38000000e+02 -4.61250000e+01 1.00207500e+20]] 

ausgibt, wird so sowohl Octave und numpy bietet die gleiche Antwort, aber es ist sehr verschieden von Matlab. Der erste Punkt ist, dass V1! = V2, was nicht richtig erscheint. Der andere Punkt ist, dass, obwohl die nicht diagonalen Elemente viele Größenordnungen kleiner sind als die diagonalen Elemente, dies in meinem Algorithmus ein Problem verursacht.

Kennt jemand, wie sich numpy und Octave auf diese Weise verhalten?

Answer 1

Sie verwenden interne Doppel, die nur etwa 15 Stellen Genauigkeit haben. Ihre mathematischen Operationen überschreiten das wahrscheinlich, was zu den mathematisch falschen Ergebnissen führt.

Verdient eine Lese: http://floating-point-gui.de/

Edit: Vom docs entnehme ich, dass es keine höhere Präzision für Numpy zur Verfügung steht. Es scheint, dass SymPy obwohl may give you the needed precision - wenn diese Bibliothek auch für Sie arbeitet.

Answer 2

Denn alles, was es wert ist, bekomme ich die gleichen Ergebnisse zu Matlab auf einem 64-Bit-System:

[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 0.00000000e+00 0.00000000e+00] 
[ -8.00000000e+00 0.00000000e+00 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]] 
[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 0.00000000e+00 0.00000000e+00] 
[ -8.00000000e+00 0.00000000e+00 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]] 
[[ 1.00207500e+20 0.00000000e+00 -8.00000000e+00 0.00000000e+00] 

Wenn Sie mit einem 32-Bit-System (oder haben Sie eine 32-Bit-Version von Python und numpy, die auf einem 64-Bit-System installiert sind), stoßen Sie auf Präzisionsprobleme und erhalten eine andere Antwort, wie von @Lucero unten angegeben. Sie könnten versuchen, in diesem Fall explizit 64-Bit-Gleitkommazahlen anzugeben (die Operationen sind jedoch langsamer). Versuchen Sie beispielsweise, np.array(..., dtype=np.float64) zu verwenden.

Wenn Sie denken, Sie zusätzliche precison benötigen, können Sie np.longdouble (Selben wie auf einem 64-Bit-System und np.float96 auf 32-Bit) verwenden, aber dies kann nicht auf allen Plattformen und viele Funktionen der linearen Algebra unterstützt wird gestutzt Dinge zurück auf die native Präzision.

Auch, welche BLAS-Bibliothek verwenden Sie? Die Ergebnisse von NUMPY und OCTAVE sind wahrscheinlich gleich, da sie dieselbe BLAS-Bibliothek verwenden.

Schließlich können Sie Ihre numpy Code unten vereinfachen:

import numpy as np 
A = np.array([[1,  -0.015, -0.025, -0.035], 
       [0.015, 1,  0.035, -0.025], 
       [0.025, -0.035, 1,  0.015], 
       [0.035, 0.025, -0.015, 1]]) 
P = np.eye(4)*1e20 
print A.dot(P.dot(A.T)) 
print A.dot(P).dot(A.T) 

Answer 3

np.einsum viel näher an den MATLAB ist

In [1299]: print(np.einsum('ij,jk,lk',A,P,A)) 
[[ 1.00207500e+20 0.00000000e+00 -5.07812500e-02 0.00000000e+00] 
[ 0.00000000e+00 1.00207500e+20 5.46875000e-02 0.00000000e+00] 
[ -5.46875000e-02 5.46875000e-02 1.00207500e+20 0.00000000e+00] 
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00207500e+20]] 

Die aus diagonalen Terme in Zeile und Spalte 2 sind unterschiedlich, aber Es hat die gleichen 0s woanders.

Mit dem Doppelpunkt erzeugt P.dot(A.T) Rundungsfehler, wenn es die Produkte hinzufügt. Diese propagieren in die nächste dot. einsum generiert alle Produkte mit nur einer Summierung. Ich vermute, dass der MATLAB-Interpreter diese Situation erkennt und eine spezielle Berechnung durchführt, um die Rundungsfehler zu minimieren.

Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix - ist eine neuere Frage mit der gleichen Erklärung.